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We develop an immersed boundary (IB) method for modeling flows around fixed or mov¬ 
ing rigid bodies that is suitable for a broad range of Reynolds numbers, including steady 
Stokes flow. The spatio-temporal discretization of the fluid equations is based on a stan¬ 
dard staggered-grid approach. Fluid-body interaction is handled usiirg Peskin’s IB method; 
however, unlike existing IB approaches to such problems, we do not rely on penalty or 
fractional-step formulatioirs. Instead, we use an unsplit scheme that ensures the no-slip 
constraint is enforced exactly in terms of the Lagrangian velocity held evaluated at the IB 
markers. Fractional-step approaches, by contrast, can impose such constraints only approx¬ 
imately, which can lead to penetration of the flow into the body, and are inconsistent for 
steady Stokes flow. Imposing no-slip constraints exactly requires the solution of a large 
linear system that includes the fluid velocity and pressure as well as Lagrange multiplier 
forces that impose the motion of the body. The principal contribution of this paper is that 
it develops an efficient preconditioner for this exactly constrained IB formulation which is 
based on an analytical approximation to the Schur complement. This approach is enabled 
by the near translational and rotational invariance of Peskin’s IB method. We demonstrate 
that only a few cycles of a geometric multigrid method for the fluid equations are required 
in each application of the preconditioner, and we demonstrate robust convergence of the 
overall Krylov solver despite the approximations made in the preconditioner. We empiri¬ 
cally observe that to control the condition number of the coupled linear system while also 
keeping the rigid structure impermeable to fluid, we need to place the immersed boundary 
markers at a distance of about two grid spacings, which is significantly larger from what 
has been recommended in the literature for elastic bodies. We demonstrate the advantage 
of our monolithic solver over split solvers by computing the steady state flow through a 
two-dimensional nozzle at several Reynolds numbers. We apply the method to a number of 
benchmark problems at zero and finite Reynolds numbers, and we demonstrate first-order 
convergence of the method to several analytical solutions and benchmark computations. 
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I. INTRODUCTION 

A large number of numerical methods have been developed to simulate interactions between 
fluid flows and immersed bodies. For rigid bodies or bodies with prescribed kinematics, many 
of these approaches [IHS] are based on the immersed boundary (IB) method of Peskin [6]. The 
simplicity, flexibility, and power of the IB method for handling a broad range of fluid-structure 
interaction problems was demonstrated by Bhalla et al. [2]. In that study, the authors showed 
that the IB method can be used to model complex flows around rigid bodies moving with specified 
kinematics (e.g., swimming fish or beating flagella) as well as to compute the motion of freely 
moving bodies driven by flow. In the approach of Bhalla et ai, as well as those of others [DEHS], 
the rigidity constraint enforcing that the fluid follows the motions of the rigid bodies is imposed 
only approximately. Here and throughout this manuscript, when we refer to the no slip condition, 
we mean the requirement that the interpolated fluid velocity exactly match the rigid body velocity 
at the positions of the IB marker points. In this work, we develop an effective solution approach 
to an IB formulation of this problem that exactly enforces both the incompressibility and no-slip 
constraints, thus substantially improving upon a large number of existing techniques. 

A simple approach to implementing rigid bodies using the traditional IB method is to use 
stiff springs to attach markers that discretize the body to tether points constrained to move as a 
rigid body [7]. This penalty-spring approach leads to numerical stiffness and, when the forces are 
handled explicitly, requires very small time steps. For this reason, a number of direct forcing IB 
methods [8] have been developed that aim to constrain the flow inside the rigid body by treating the 
fluid-body force as a Lagrange multiplier A enforcing a no-slip constraint at the locations of the IB 
markers. However, to our knowledge, all existing direct forcing IB methods use some form of time 
step splitting to separate the coupled fluid-body problem into more manageable pieces. The basic 
idea behind these approaches is first to solve a simpler system in which a number of the constraints 
(e.g., incompressibility, or no-slip along the fluid-body interface) are ignored. The solution of the 
unconstrained problem is then projected onto the constraints, which yields estimates of the true 
Lagrange multipliers. In most existing methods, the fluid solver uses a fractional time stepping 
scheme, such as a version of Chorin’s projection method, to separate the velocity update from the 
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pressure update PEHS]. Taira and Colonius also use a fractional time-stepping approach in which 
they split the velocity from the Lagrange multipliers (tt, A). They obtain approximations to (vr, A) 
in a manner similar to that in a standard projection method for the incompressible Navier-Stokes 
equations. A modified Poisson-type problem (see (26) in [3]) determines the Lagrange multipliers 
and is solved using an unpreconditioned conjugate gradient method. The method developed in Ref. 
[ 2 ] avoids the pressure-velocity splitting and instead uses a combined iterative Stokes solver, and 
in Ref. [3j (see supplementary material), periodic boundary conditions are applied, which allows 
for the use of a pseudo-spectral method. In both works, however, time step splitting is still used to 
separate the computation of the rigidity constraint forces from the updates to the fluid variables. 
In the approach described in the supplementary material to Ref. [ 3 ], the projection step of the 
solution onto the rigidity constraint is performed twice in a predictor-corrector framework, which 
improves the imposition of the constraint; however, this approach does not control the accuracy of 
the approximation of the constraint forces. Curet et al. [9] and Ardekani et al. m go a step closer 
in the direction of exactly enforcing the rigidity constraint by iterating the correction until the 
relative slip between the desired and imposed kinematics inside the rigid body reaches a relatively 
loose tolerance of 1%. The scheme used in Ref. [9] is essentially a fixed-point (Richardson) 
iteration for the constrained fluid problem, which uses splitting to separate the update of the 
Lagrange multipliers from a fluid update based on the SIMPLER scheme [m. Unlike the approach 
developed here, fixed point iterations based on splitting are not guaranteed to converge, yet alone 
converge rapidly, especially in the steady Stokes regime for tight solver tolerances. 

An alternative view of direct forcing methods that use time step splitting is that they are penalty 
methods for the unsplit problem, in which the penalty parameter is related to the time step size. 
Such approaches inherently rely on inertia and implicitly assume that fluid velocity has memory. 
Consequently, all such splitting methods fail in the steady Stokes limit. Furthermore, even at 
finite Reynolds numbers, methods based on splitting cannot exactly satisfy the no-slip constraint 
at fluid-body interfaces. Such methods can thereby produce undesirable artifacts in the solution, 
such as penetration of the flow through a rigid obstacle. It is therefore desirable to develop a 
numerical method that solves for velocity, pressure, and fluid-body forces in a single step with 
controlled accuracy and reasonable computational complexity. 

The goal of this work is to develop an effective IB method for rigid bodies that does not rely on 
any splitting. Our method is thus applicable over a broad range of Reynolds numbers, including 
steady Stokes flow, and is able to impose rigidity constraints exactly. This approach requires us 
to solve large linear systems for velocity, pressure, and fluid-body interaction forces. This linear 
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system is not new. For example, (13) in Ref. [3] is essentially the same system of equations that we 
study here. The primary contributions of this work are that we do not rely on any approximations 
when solving this linear system, and that we develop an effective preconditioner based on an 
approximation of the Schur complement that allows us to solve Q. The resulting method has a 
computational complexity that is only a few times larger than the corresponding problem in the 
absence of rigid bodies. In the context of steady Stokes flows, a rigid-body IB formulation very 
similar to the one we use here has been developed by Bringley and Peskin [12]; however, that 
formulation relies on periodic boundary conditions, and uses a very different spatial discretization 
and solution methodology from the approach we describe here. Our approach can readily handle a 
broad range of specified boundary conditions. In both Refs. [12] and a very recent work by Stein 
et al. on a higher-order IB smooth extension method for scalar (e.g., Poisson) equations [8], the 
Schur complement is formed densely in an expensive pre-computation stage. By contrast, in the 
method proposed here we build a simple physics-based approximation of the Schur complement 
that can be computed “on the fly” in a scalable and efficient manner. 

Our basic solution approach is to use a preconditioned Krylov solver for the fully constrained 
fluid problem, as has been done for some time in the context of finite element methods for fluid 
flows interacting with elastic bodies [nun]. A key difficulty that we address in this work is the 
development of an efficient preconditioner for the constrained formulation. To do so, we construct 
an analytical approximation of the Schur complement (i.e., the mobility matrix) corresponding to 
Lagrangian rigidity forces (i.e., Lagrange multipliers) enforcing the no-slip condition at the positions 
of the IB markers. We rely on the near translational and rotational invariance of Peskin’s IB method 
to approximate the Schur complement, following techniques commonly used for suspensions of 
rigid spheres in steady Stokes flow such as Stokesian dynamics [13 US], bead methods for rigid 
macromolecules [iiHnn] and the method of regularized Stokeslets [2TII23] . In fact, as we explain 
herein, many of the techniques developed in the context of steady Stokes flow can be used with the 
IB method both at zero and also, perhaps more surprisingly, finite Reynolds numbers. 

The method we develop offers an attractive alternative to existing techniques in the context of 
steady or nearly-steady Stokes flow of suspensions of rigid particles. To our knowledge, most other 
approaches tailored to the steady Stokes limit rely on Green’s functions for Stokes flow to eliminate 
the (Eulerian) fluid degrees of freedom and solve only for the (Lagrangian) degrees of freedom 
associated to the surface of the body. Because these approaches rely on the availability of analytical 
solutions, handling non-trivial boundary conditions (e.g., bounded systems) is complicated [24] and 
has to be done on a case-by-case basis [25H3T] . By contrast, in the method developed here, analytical 
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Green’s functions are replaced by an “on the fly” computation that may be carried out by a standard 
finite-volume, finite-difference, or finite-element fluid solver Such solvers can readily handle 
nontrivial boundary conditions. Furthermore, suspensions at small but nonzero Reynolds numbers 
can be handled without any extra work. Additionally, we avoid uncontrolled approximations relying 
on truncations of multipole expansions to a fixed order [I51I31H36], and we can seamlessly handle 
arbitrary body shapes and deformation kinematics. For problems involving active m particles, it is 
straightforward to add osmo- or electro-phoretic coupling between the fluid flow and additional fluid 
variables such as the electric potential or the concentration of charged ions or chemical reactants. 
Lastly, in the spirit of fluctuating hydrodynamics |38H4U] . it is straightforward to generate the 
stochastic increments required to simulate the Brownian motion of small rigid particles suspended 
in a fluid by including a fluctuating stress in the fluid equations. We also point out that our method 
also has some disadvantages compared to methods such as boundary integral or boundary element 
methods. Notably, it requires filling the domain with a dense uniform fluid grid, which is expensive 
at low densities. It is also a low-order method that cannot compute solutions as accurately as 
spectral boundary integral formulations. We do believe, nevertheless, that the method developed 
here offers a good compromise between accuracy, efficiency, scalabilty, flexibility and extensibility, 
compared to other more specialized formulations. 

II. SEMI-CONTINUUM FORMULATION 

Our notation uses the following conventions where possible. Vectors (including multi-vectors), 
matrices, and operators are bolded, but when fully indexed down to a scalar quantity we no longer 
bold the symbol; matrices and operators are also scripted. We denote Eulerian quantities with 
lowercase letters, and the corresponding Lagrangian quantity with the same capital letter. We use 
the Latin indexes i,j, k, I, m to denote a specific fluid grid point or IB marker (i.e., physical location 
with which degrees of freedom are associated), the indices p, q, r, s, t to denote a specific body in the 
multibody context, and Greek superscripts a ,/?,7 to denote specihc Cartesian components. For 
example, v denotes fluid velocity (either continuum or discrete), with being the fluid velocity 
in direction a associated with the face center k, and V denotes the velocity of all IB markers, 
with being the velocity of marker i along direction a. Our formulation is easily extended to a 
collection of rigid bodies, but for simplicity of presentation, we focus on the case of a single body. 

^ In this work, we use a staggered-grid discretization on a uniform grid combined with multigrid-preconditioned 
Stokes solvers |321133 |. 
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We consider a region 2? C (d = 2 or 3) that contains a single rigid body Q (ZV immersed in 
a fluid of density p and shear viscosity p. The computational domain D could be a periodic region 
(topological torus), a finite box, an infinite domain, or some combination thereof, and we will 
implicitly assume that some consistent set of boundary conditions are prescribed on its boundary 
dD even though we will not explicitly write this in the formulation. We require that the linear 
velocity of a given reference point (e.g., the center of mass of the body) U{t) and the angular 
velocity fl{t) of the body are specified functions of time, and without loss of generality, we assume 
that the rigid body is at rest In addition to features of the fluid flow, typical quantities of 
interest are the total drag force F{t) and total drag torque T(t) between the fluid and the body. 
Another closely related problem to which IB methods can be extended is the case when the motion 
of the rigid body (i.e., U{t) and r2(t)) is not known but the body is subject to specified external 
force F{t) and torque T{t). For example, in the sedimentation of rigid particles in suspension, the 
external force is gravity and the external torque is zero. Handling this free kinematics problem 
mm requires a nontrivial extension of our formulation and numerical algorithm. 

In the immersed boundary (IB) method |6l l42( H3]. the velocity field v(r, t) is extended over the 
whole domain F, including the body interior. The body is discretized using a collection of markers, 
which is a set of N points that cover the interior of the body and at which the interaction between 
the body and the fluid is localized. For example, the markers could be the nodes of a triangular 
{d = 2) or tetrahedral {d = 3) mesh used to discretize an illustration of such a volume grid of 
markers discretizing a rigid disk immersed in steady Stokes flow is shown in the left panel of Fig. 

In the case of Stokes flow, the specihcation of a no slip condition on the boundary of a rigid body 
is sufficient to ensure rigidity of the fluid inside the body [H]. Therefore, for Stokes flow, the grid 
of markers does not need to extend over the volume of the body and can instead be limited to the 
surface of the rigid body, thus substantially reducing the number of markers required to represent 
the body. In this case, the markers could be the nodes of a triangulation {d = 3) of the surface of 
the body; an illnstration of such a surface grid of markers is shown in the right panel of Fig. 

We discuss the differences between a volume and a surface grid of markers in Section [VII 

The traditional IB method is concerned with the motion of elastic (flexible) bodies in fluid 
flow, and the collection of markers can be viewed as a set of quadrature points used to discretize 
integrals over the moving body. The elastic body forces are most easily computed in a Lagrangian 

^ The case of more general specified kinematics is a straightforward generalization and does not incur any additional 
mathematical or algorithmic complexity [S] [41]. 
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coordinate system attached to the deforming body, and the relative positions of the markers in the 
fixed Eulerian frame of reference generally vary in time. For a rigid body, however, the relative 
positions of the markers do not change, and it is not necessary to introduce two distinct coordinate 
frames. Instead, we use the same Cartesian coordinate system to describe points in the fluid domain 
and in the body; the positions of the N markers in this fixed frame of reference will be denoted 
with R = {/?!,..., Rn}, where i? C for volume meshes or R C dQ for surface meshes. 



Figure 1: Two-dimensional Steady Stokes flow past a periodic column of circular cylinders (disks) at zero 
Reynolds number obtained using our rigid-body IB method (the same setup is also studied at finite Reynolds 


numbers in Section VIIF). The markers used to mediate the fluid-body interaction are shown as small 
colored circles. The Lagrangian constraint forces A that keep the markers at their fixed locations are shown 
as colored vectors; the color of the vectors and the corresponding marker i are based on the magnitude of 
the constraint force (see color bar). The fluid velocity field is shown as a vector field (black arrows) in the 
vicinity and the interior of the body; further from the body, flow streamlines are shown as solid blue lines. 
The magnitude of the Eulerian constraint force «SA is shown as a gray color plot (see greyscale bar). (Left 
panel) A volume marker grid of 121 markers is used to discretize the disk. The majority of the constraint 
forces are seen to act near the surface of the body, but nontrivial constraint forces are seen also in the 
interior of the body. (Right panel) A surface grid of 39 markers is used to discretize the disk, which strictly 
localizes the constraint forces to the surface of the body. 


In the standard IB method for flexible immersed bodies, elastic forces are computed in the 
Lagrangian frame and then spread to the fluid in the neighborhood of the markers using a regu¬ 
larized delta function 6a (r) that integrates to unity and converges to a Dirac delta function as the 
regularization width a —)• 0. The regularization length scale a is typically chosen to be on the order 
of the spacing between the markers (as well as the lattice spacing of the grid used to discretize the 
fluid equations), as we discuss in more detail later. In turn, the motion of the markers is specihed 
to follow the velocity of the fluid interpolated at the positions of the markers. 

The key difference between an elastic and a rigid body is that, for a rigid object, the motion of 
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the markers is known (e.g, they are fixed in place or move with a specified velocity) and the body 
forces are unknown and must be determined within each time step. To obtain the fluid-marker 
interaction forces A (t) = {Aj (t),..., Ajv (t)} that constrain the motion of the N markers, we 
solve for the Eulerian velocity field v (r, t), the Eulerian pressure field tt (r, t), and the Lagrangian 
constraint forces Aj (t) the system 

N 

p {dtv + V • Vv) + Vtt = p'V'^v + ^ Ai6a {Ri - r ), 

i=l 

V-u = 0, 

Vi = j da{Ri-r)v{r,t) dr = 0, i = (1) 

along with suitable boundary conditions. In the case of steady Stokes flow, we set p = 0. The first 
two equations are the incompressible Navier-Stokes equations with an Eulerian constraint force 

N 

A (r, i) = ^ Aida {Ri - r). 

i=l 

The last condition is the rigidity constraint that requires that the Eulerian velocity averaged around 
the position of marker i must match the known marker velocity Vi. This constraint enforces a 
regularized no-slip condition at the locations of the IB markers, which is a numerical approximation 
of the true no-slip condition on the surface (or interior) of the body. Observe that flow may still 
penetrate the body in-between the markers and this leads to a well-known small but nonzero “leak” 
in the traditional Peskin IB method. This leak can be greatly reduced by adopting a staggered- 
grid formulation [44j . as done in the present work. Other more specialized approaches to reducing 
spurious fluxes in the IB method have been developed [iHIHT] . but will not be considered in this 
work. 

Notice that for zero Reynolds number, the semi-continuum formulation Q is closely related 
to the popular method of regularized Stokeslets, which solves a similar system of equations for A 
HUES]. The key difference ^ is that in the method of regularized Stokeslets, the fluid equations 
are eliminated using analytic Green’s functions; this necessitates that nontrivial pre-computations 
of these Green’s functions be performed for each type of boundary condition [301 EH- 

® Another important difference is that we follow Peskin and use the regularized delta function both for spreading 
and interpolation (this ensures energy conservation in the formulation [^), whereas in the method of Regularized 
Stokeslets only the spreading uses a regularized delta function. Our choice ensures that the linear system we solve 
is symmetric and positive semi-definite, which is crucial if one wishes to account for Brownian motion and thermal 
fluctuations. 
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In this work, we treat 0 as the primary continuum formulation of the problem. This is a semi¬ 
continuum formulation in which the rigid body is represented as a discrete collection of markers but 
the fluid description is kept as a continuum, which implies that different discretizations of the fluid 
equations are possible. One can, in principle, try to write a fully continuum formulation in which 
the discrete set of rigidity forces A are replaced by a continuum force density field \{R G 0,t). 
The well-posedness and stability of such a fully continuum formulation is mathematically delicate, 
however, and there can be subtle differences between weak and strong interpretations of the equa¬ 
tions. To appreciate this, observe that if each component of the velocity is discretized with Nf 
degrees of freedom, it cannot in general be possible to constrain the velocity strongly at more than 
Nf points (markers). By contrast, in our strong formulation Q, the velocity is infinite dimensional 
but it is only constrained in the vicinity of a finite number of markers. Therefore, the problem 
0 is always well posed and is directly amenable to numerical discretization and solution, at least 
when it is well-conditioned. As we show in this work, the conditioning of the fully discrete problem 
is controlled by the relationship between the regularization length a and the marker spacing. 

The physical interpretation of the constraint forces Aj depends on details of the marker grid 
and the type of the problem under consideration. For fully continuum formulations, in which the 
fluid-body interaction is represented solely as a surface force density, the force A* can be interpreted 
as the integral of the traction (normal component of the fluid stress tensor) over a surface area 
associated with marker i. Such a formulation is appropriate, for example, for steady Stokes flow. In 
particular, for steady Stokes flow our method can be seen as a discretized and regularized first-kind 
integral formulation in which Green’s functions are computed by the fluid solver. This approach 
is different from the method of regularized Stokeslets, in which regularized Green’s functions must 
be computed analytically 

For cases in which markers are placed on both the snrface and the interior of a rigid body, 
the precise physical interpretation of the volume force density, and thus of A, is delicate even for 
steady Stokes flow. Notably, observe that the splitting between a volume constraint force density 
and the gradient of the pressnre is not unique because the pressure inside a rigid body cannot be 
determined uniquely. Specihcally, only the component of the constraint force density projected 
onto the space of divergence-free vector fields is nniqnely determined. In the presence of finite 
inertia and a density mismatch between the fluid and the moving rigid bodies, the inertial terms 
in 0 need to be modified in the interior of the body |3S]- Furthermore, sufficiently many markers 
in the interior of the body are required to prevent spnrious angular momentum being generated 
by motions of the fluid inside the body [T]. We do not discuss these physical issues in this work 
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because they do not affect the numerical algorithm, and because we restrict our numerical studies 
to flow past stationary rigid bodies, for which the fluid-body interaction force is localized to the 
surface of the body in the continuum limit. 

III. DISCRETE FORMULATION 

The spatial discretization of the fluid equation uses a uniform Cartesian grid with grid spacing 
h and is based on a second-order accurate staggered-grid finite-difference discretization, in which 
vector-valued quantities, including velocities and forces, are represented on the faces of the Carte¬ 
sian grid cells, and scalar-valued quantities, including the pressure, are represented at the centers 
of the grid cells [2113511121113]. Our implicit-explicit temporal discretization of the Navier-Stokes 
equation is standard and summarized in prior work; see for example the work of Griffith |43) . The 
key features are that we treat advection explicitly using a predictor-corrector approach, and that 
we treat viscosity implicitly, using either the backward Euler or the implicit midpoint method. For 
steady Stokes flow, no temporal discretization required, although one can also think of this case as 
corresponding to a backward Euler discretization of the time-dependent problem with a very large 
time step size At. A key dimensionless quantity is the viscous CEL number (3 = uAt/h'^, where 
the kinematic viscosity is v = rj/p. If /3 is small, the pressure and velocity are weakly coupled, but 
for large (3, and in particular for the steady Stokes limit /3 —>■ oo, the coupling between the velocity 
and pressure equations is strong. 

We do not use a fractional time-stepping scheme (i.e., a projection method) to split the pressure 
and velocity updates; instead, the pressure is treated as a Lagrange multiplier that enforces the 
incompressibility and must be determined together with the velocity at the end of the time step 
|32| : except in special cases, this is necessary for small Reynolds number flows. This approach also 
greatly aids with imposing stress boundary conditions [32]. The constraint force X{r,t) is treated 
analogously to the pressure, i.e., as a Lagrange multiplier. Whereas the role of the pressure is to 
enforce the incompressibility constraint, A enforces the rigidity constraint. Like the pressure, A is 
an unknown that must be solved for in this formulation. 

A. Force spreading and velocity interpolation 

In the fully discrete formulation of the fluid-body coupling, we replace spatial integrals by sums 
over fluid or body grid points in the semi-continuum formulation Q . The regularized delta function 
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is discretized using a tensor product in d-dimensional space (see [35] for more details), 

d 

5a{r) = h~‘^ 0a (ra) , 
a=l 

where h'^ is the volume of a grid cell. The one-dimensional kernel function (pa is chosen based on 
numerical considerations of efficiency and maximized approximate translational invariance |6|. In 


this work, for reasons that will become clear in Section IV, we prefer to use a kernel that maximizes 
translational and rotational invariance (i.e., improves grid-invariance). We therefore use the smooth 
(three-times differentiable) six-point kernel recently described by Bao et al. |48) . This kernel is more 
expensive than the traditional four-point kernel [6] because it increases the support of the kernel to 
6 ^ = 36 grid points in two dimensions and 6^ = 216 grid points in three dimensions; however, this 
cost is justified because the new six-point kernel improves the translational invariance by orders of 
magnitude compared to other standard IB kernel functions |48j . 

The interaction between the fluid and the rigid body is mediated through two crucial operations. 
The discrete velocity-interpolation operator PT averages velocities on the staggered grid in the 
neighborhood of marker i via 


{Jv)^ = Y,vt <PaiRi-rt), 

k 

where the sum is taken over faces k of the grid, a indexes coordinate directions (x, y, z) as a 
superscript, and is the position of the center of the grid face k in the direction a. The discrete 
force-spreading operator S spreads forces from the markers to the faces of the staggered grid via 


{S\)t = h-‘‘Y.K4>a{Ri-rt), ( 2 ) 

i 

where now the sum is over the markers that define the configuration of the rigid body. These 
operators are adjoint with respect to a suitably-defined inner product, PT = S* = , which 

ensures conservation of energy (Hj. Extensions of the basic interpolation and spreading operators 
to account for the presence of physical boundary conditions are described in Appendix [D| 


B. Rigidly-constrained Stokes problem 

At every stage of the temporal integrator, we need to solve a linear system of the form 


A g -s 


V 


9 

o 

o 

1 


TT 

= 

h = 0 

-J" 0 0 


A 


W = 0 


( 3 ) 
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which is the focus of this work. The right-hand side g includes all remaining fluid forcing terms, 
explicit contributions from previous time steps or stages, boundary conditions, etc. Here, Q is 
the discrete gradient operator, T> = is the discrete divergence operator, and A. is the vector 
equivalent of the familiar screened Poisson (or Helmholtz) operator 


A=^X-^C 

At 


with K = 1 for the backward Euler method or for steady Stokes, and k = 1/2 for the implicit 
midpoint rule. Here is the dimensionless vector Laplacian operator, which takes into account 
boundary conditions for velocity such as no-slip boundaries. Since the viscosity appears multiplied 
by the coefficient k, we will henceforth absorb this coefficient into the viscosity, rj •(— ng, which 
allows us to assume, without loss of generality, that k=\ and to write the fluid operator in the 
form 


A = gh-^ - £,) . (4) 

We remark that making the (3,3) block in the matrix in ([^ non-zero (i.e., regularizing the 
saddle-point system) is closely related to solving the Brinkman equations [lU] for flow through 
a permeable or porous body suspended in fluid [1]. In particular, by making the (3,3) block a 
diagonal matrix with suitable diagonal elements, one can consistently discretize the Brinkman 
equations. Such regularization greatly simplifies the numerical linear algebra except, of course, 
when the permeability of the body is so small that it effectively acts as an impermeable body. In 
this work, we focus on developing a solver for @ that is effective even when there is no regularization 
(permeability), and even when the matrix A is the discretization of an elliptic operator, as is the 
case in the steady Stokes regime. This is the hardest case to consider, and a solver that is robust 
in this case will be able to handle the easier cases of finite Reynolds number or permeable bodies 
with ease. 

It is worth noticing the structure of the linear system ([^. First, observe that the system is 
symmetric, at least if only simple boundary conditions such as periodic or no-slip boundaries are 
present |32| . In the top 1x1 block, ^ 0 is a symmetric positive-semidefinite (SPD) matrix. The 

top left 2x2 block represents the familiar saddle-point problem arising when solving the Navier- 
Stokes or Stokes equations in the absence of a rigid body [32] . The whole system is a saddle-point 
problem for the fluid variables and for A, in which the top-left block is the Stokes saddle-point 
matrix. 
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C. Mobility matrix 


We can formally solve ^ through a Schur complement approach, as described in more detail in 
Section |Vj For increased generality, which will be useful when discussing preconditioners, we allow 
the right hand side to be general and, in particular, do not assume that h and W are zero. 

First, we solve the unconstrained fluid equation for pressure and velocity 


A 

G 


V 


SA + g 

-T> 

0 


TT 


h 


where we recall that A. = rjh ^ (/? . The solution can be written as u = X ^ (<SA + g) + 

Cp^h, where C~^ is the standard Stokes solution operator for divergence-free flow {h = 0), given 
by 


X-^ = A-^ - A {T>A-^G) ^ T>A \ 


( 6 ) 


where we have assumed for now that ^ ^ is invertible. For a periodic system, the discrete operators 
commute, and we can write 


C-^ =TA-^ = (l-giVG)-^v'^A-\ (7) 

where 7^ is the Helmholtz projection onto the space of divergence-free vector fields. We never 
explicitly compute or form X~^; rather, we solve the Stokes velocity-pressure subsystems using the 
projection-method based preconditioner developed by Griffith [32]. Let us define v = C~^ f+Cp^h 
to be the solution of the unconstrained Stokes 

A G 
-T> 0 

giving V = V + C^^SA. 

Next, we plug the velocity v into the rigidity constraint, J^v = — W, to obtain 


problem 


V 


1- 

1_ 

TT 


h 


( 8 ) 


MA = -{W + Jv ), 


(9) 


where the Schur complement or marker mobility matrix is 

M = JC = S*C~^S. (10) 

The mobility matrix A4 ^ 0 is SPD and has dimensions dN x dN, and the dx d block Ad-ij relates 
the force applied at marker j to the velocity induced at marker i. Our approach to obtaining an 
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efficient algorithm for the constraind fluid-solid system is to develop a method for approximating 
the marker mobility matrix A4 in a simple and efficient way that leads to robust preconditioners 
for solving the mobility subproblem Q; see Section 

Observe that the conditioning of the saddle-point system is controlled by the conditioning 
of A4. In particular, if the (non-negative) eigenvalues of A4 are bounded away from zero, then 
there will be a unique solution to the saddle-point system. If this bound is uniform as the grid 
is refined, then the problem is well-posed and will satisfy a stability criterion similar to the well- 
known Ladyzenskaja-Babuska-Brezzi (LBB) condition for the Stokes saddle-point problem Q. We 
investigate the spectrum of the the marker mobility matrix numerically in Section |Vj In practice, 
there may be some nearly zero eigenvalues of the matrix A4 corresponding to physical (rather than 
numerical) null modes. An example is a sphere discretized with markers on the surface: we know 
that a uniform compression of the sphere will not cause any effect because of the incompressibility 
of the fluid hlling the sphere. This compression mode corresponds to a null-vector for the constraint 
forces A; it poses no difficulties in principle because the right-hand side in Q is always in the range 
of A4. Of course, when a discrete set of markers is placed on the sphere, the rotational symmetry 
will be broken and the corresponding mode will have a small but nonzero eigenvalue, which can 
lead to numerical difficulties if not handled with care. 

D. Periodic steady Stokes flow 

In the time-dependent context, /3 is finite, and it is easy to see that ^ ^ 0 is invertible. The 
same happens even for steady Stokes flow if at least one of the boundaries is a no-slip boundary. 
In the case of periodic steady Stokes flow, however, A. = —rih~‘^Cy has in its range vectors that 
sum to zero, because no nonzero total force can be applied on a periodic domain. This means that 
a solvability condition is 

N 

{SK + g) = vol”^ ^ Aj -h {g) = voffi^'^A -F {g) = 0, 

i=l 

where () denotes an average over the whole system, 1 is a vector of ones, and vol is the volume 
of the domain. This is an additional constraint that must be added to the constrained Stokes 
system Q for a periodic domain the steady Stokes case. In this approach, the solution has an 
indeterminate mean velocity {v) because momentum is not conserved. This sort of approach is 
followed for a scalar (reaction-diffusion) equivalent of Q in the Appendix of Ref. m, for the 
traditional Beskin IB method in Ref. [7], and for a higher-order IB method in [8]. 


IV 
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Here, we instead impose the mean velocity {v) = 0 and ensure that the total force applied to 
the fluid sums to zero, i.e., we enforce momentum conservation. Specifically, for the special case of 
periodic steady Stokes, we solve the system 


A Q 

- {S - yoI~A'^) 


V 


9 

o 

1 

0 


TT 

= 

h 

-J 0 

0 


A 


W 


together with the constraint {v) = 0, where we assume that (g) = 0 for consistency. This change 
amounts to simply redefining the spreading operator to subtract the total applied force on the 
markers as a uniform force density, S S — vol~^l^. This can be justified by considering the 
unit cell to be part of an infinite periodic system in which there is an externally applied constant 
pressure gradient, which is balanced by the drag forces on the bodies so as to ensure that the 
domain as a whole is in force balance [SIHS3]- 

IV. APPROXIMATING THE MOBILITY MATRIX 

A key element in the preconditioned Krylov solver for Q that we describe in Section |V] is an 
approximate solver for the mobility subproblem Q. The success of this approximate solver, i.e., 
the accuracy with which we can approximate the Schur complement of the saddle-point problem 
(§, is crucial to an effective linear solver and one of the key contributions of this work. 

Because it involves the inverse Stokes operator C~^, the actual Schur complement A4 = 
S*C~^S cannot be formed efficiently. Instead of forming the true mobility matrix, we instead 
approximate A1 « A4 by a dense but low-rank approximate mobility matrix A4 given by simple 
analytical approximations. To achieve this, we use two key ideas: 

1. We ignore the specifics of the boundary conditions and assume that the structure is immersed 
in an inhnite domain at rest at infinity (in three dimensions) or in a hnite periodic domain 
(in two dimensions). This implies that the Krylov solver for ([^ must handle the boundary 
conditions. 

2. We assume that the IB spatial discretization is translationally and rotationally invariant; 
that is, A4 does not depend on the exact position and orientation of the body relative to the 
underlying fluid grid. This implies that the Krylov solver must handle any grid-dependence 
in the solution. 
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The first idea, to ignore the boundary conditions in the preconditioner, has worked well in the 
context of solving the Stokes system Q. Namely, a simple but effective approximation of the 
inverse of the Schur complement for (j^, can be constructed by assuming that the 

domain is periodic so that the finite difference operators commute, and thus the Schur complement 
degenerates to a diagonal or nearly-diagonal mass matrix [sa Ea [51]. The second idea, to make 
use of the near grid invariance of Peskin’s regularized kernel functions, has previously been used 
successfully in implicit immersed-boundary methods by Ceniceros et al. [55]. Note that for certain 
choices of the kernel function, the assumption of grid invariance can be a very good approximation 
to reality; here, we rely on the recently-developed six-point kernel [48j . which has excellent grid 
invariance and relatively compact support. 

In the remainder of this section, we explain how we compute the entries in A4 in three dimen¬ 
sions, assuming an unbounded fluid at rest at infinity. The details for two dimensions are given in 
Appendix [B] and are similar in nature, except for complications for two-dimensional steady Stokes 
flow resulting from the well-known Stokes paradox. 

The mobility matrix Ad is a symmetric block matrix built from N x N blocks of size d x d. 
The block Adjj corresponding to markers i and j relates a force applied at marker j to the velocity 
induced at marker i. Our basic assumption is that Aiij does not depend on the actual position 
of the markers relative to the fluid grid, but rather only depends on the distance between the two 
markers and on the viscous CFL number /3 in the form 

■^ij = //3 irij)^ + 913 (rij) Tij (g) Tij, (12) 


where Vij = Hi — Rj and is the distance between the two markers, and hat denotes a unit 
vector. The functions of distance //^(r) and 5/3 (r) depend on the specific kernel chosen, the specific 
discretization of the fluid equations (in our case the staggered-grid scheme), and the the viscous 
CFL number (5. To obtain a specific form for these two functions, we empirically fit numerical 
data with functions with the proper asymptotic behavior at short and large distances between the 
markers. For this purpose, we first discuss the asymptotic properties of //^(r) and 5 / 3 (r) from a 
physical perspective. 

It is important to note that the true mobility matrix A4 is guaranteed to be SPD because of its 
structure and the adjointness of the spreading and interpolation operators. This can be ensured 
for the approximation A4 by placing positivity constraints on suitable linear combinations of the 


Fourier transforms of (r) and gjsir), which ensure that the kernel A4 (r,;. rj) given by (12) is 


SPD in the sense of integral operators. It is, however, very difficult to place such constraints on 
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empirical fits in practice, and in this work, we do not attempt to ensure Ad is SPD for all marker 
configurations. 


A. Physical Constraints 

Let us temporarily focus on the semi-continuum formulation 0 and ignore Eulerian discretiza¬ 
tion artifacts. The pairwise mobility between markers i and j for a continuum fluid is 


A4^j = J 5a{R^ - r")G{r",r')5a{R, - r') dr"dr', 

where G{r,r') is the the Green’s function for the fluid equation, i.e., v{r) = 
f G(r,r')f(r')dr', where 

Vvr - 7]V^v = f, 

V-u = 0. 


(13) 


r = 


(14) 


It is well-known that G has the same form as (12), 


G{Ri,R 2) = f {ri2)T + g{ri2)ri2Gri2. 

For steady Stokes flow (/3 —>■ oo), G = O is the well-known Oseen tensor or Stokeslet and 
corresponds to /s(r) = gsir) ~ (Svrr/r)”^. For inviscid flow, /? = 0, and we have that A. = {p/At)X 
and 0 applies, and therefore C ^ = {At/p)V is a multiple of the projection operator. For hnite 
nonzero values of (3, we can obtain G from the solution of the screened Stokes (i.e., Brinkman) 
equations (14) [23lll9l[56], and corresponds to the “Brinkmanlet” 


Mr) = ^((-) + ( 15 ) 

47rT/r \\arj ar I 

e-"*' 3 \ 3 

Airpr y \ar J ar~^ j ^ dvrrya^r^’ 

where c? = pj (pAt) = (/3/i^) ^. Note that in the steady Stokes limit, a —)• 0 and the Brinkmanlet 
becomes the Stokeslet. 


We can use (15) to construct AAij when the markers are far apart. Namely, if Vij h, then 
we may approximate the IB kernel function by a true delta function, and thus /^{r) and gp(r) are 


Observe that the regularized Stokeslet of Cortez m is similar to ( |l3[ ) but contains only one regularized delta 
function in the integrand; this makes the resulting mobility matrix asymmetric, which is unphysical. 
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well-approximated by (15). For steady Stokes flow, the interaction between markers decays like 
r~^. For finite /3, however, the viscous contribution decays exponentially fast as exp (—r/ (hy^)), 
which is consistent with the fact that markers interact via viscous forces only if they are at a 
distance not much larger than h\/~^ = \/z/Af, the typical distance that momentum diffuses during 
a time step. For nonzero Reynolds numbers, the leading order asymptotic r~^ decay of fj 3 {r) and 


g/ 3 {r) is given by the last terms on the right hand side of (15) and corresponds to the electric field 
of an electric dipole; its physical origin is in the incompressibility constraint, which instantaneously 
propagates hydrodynamic information between the markers 

For steady Stokes flow, we can say even more about the approximate form of ffs{r) and g/ 3 {r). 
As discussed in more detail by Delong et ah [38], for distances between the markers that are not 
too small compared to the regularization length a, we can approximate (13) with ( |12| ) using the 
well-known Rotne-Prager-Yamakawa (RPY) [571159] tensor for the functions // 3 (r) and 5^/3 (r). 


fRPY{r) 


gRPY{r) 


1 

Girga 

1 

bvTTya 


3a 

4r 2r^ ’ 

< 

1 _ 9r 

/ 32a ’ 

/ 

3a _ 3a® 

4r 2r® ’ 

< 

3r 

^ 32a’ 


r > 2 a, 
r < 2 a, 

r > 2 a, 
r < 2 a, 


(16) 


where a is the effective hydrodynamic radius of the specific kernel 5a, defined by (Ovra)”^ = 
f 5a(r")0(r",r')5a(r') dr”dr'. Note that for r » a the RPY tensor approaches the Oseen tensor 
and decays like r~^. A key advantage of the RPY tensor is that it guarantees that the mobility 


matrix (12) is SPD for all configurations of the markers, which is a rather nontrivial requirement 


|59j . The actual discrete pairwise mobility ATy obtained from the spatially-discrete IB method 
is well-described by the RPY tensor [38| (see Fig. [^. The only fitting parameter in the RPY 
approximation is the effective hydrodynamic radius a averaged over many positions of the marker 
relative to the underlying grid |35L 138] : for the six-point kernel used here a = 1.47h. For the 
Brinkman equation, the equivalent of the RPY tensor can be computed for r > 2a by applying 
a Faxen-like operator from the left and right on the Brinkmanlet (see Eq. (26) in Ref. jSSj); the 
resulting analytical expressions are complex and are not used in our empirical fitting. 


® In reality, of course, this information is propagated via fast sound waves and not instantaneously. 

® As summarized in Refs. |35l I38| . a ~ 1.25/i for the widely used four-point kernel [^, and a ~ 0.91/i for the 
three-point kernel | 60| . 
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B. Empirical Fits 


In this work, we use empirical fits to approximate the mobility. This is because the analytical 
approximations, such as those offered by the RPY tensor, are most appropriate for unbounded 
domains and assume the markers are far apart compared to the width of the regularized delta 
function. In numerical computations, we use a finite periodic domain, and this requires corrections 
to the analytic expressions that are difficult to model. For example, for finite (3, we find that 
the periodic corrections to the inviscid (dipole) r~^ contribution dominate over the exponentially 


decaying viscous contribution, which makes the precise form of the viscous terms in (15) irrelevant 


in practice. For r 3> h, only the asymptotically-dominant far-field terms survive, and we make an 
effort to preserve those in our fitting because the numerical results are obtained using finite systems 
and thus not reliable at large marker distances. At shorter distances, however, the discrete nature 
of the fluid solver and the IB kernel functions becomes important, and empirical fitting seems to 
be a simple yet flexible alternative to analytical computations. At the same time, we feel that is 
important to constrain the empirical fits based on known behavior at short and large distances. 

Firstly, for r h, the pairwise mobility can be well-approximated by the self-mobility (r = 0, 
corresponding to the diagonal elements Ai-a), for which we know the following facts: 

• For the steady Stokes regime (/3 —)■ oo), the diagonal elements are given by Stokes’s drag 
formula, yielding 


/oo(0) = (67rr/a) 


-1 


1/rih and 500(0) = 0, 


where we recall that a is the effective hydrodynamic radius of a marker for the particular 
spatial discretization (kernel and fluid solver). 


For the inviscid case (/? = 0), it is not hard to show that j35| 

/o(0) = ~ l3/rih and 50 ( 0 ) = 0, 

a p 


(17) 


where d = 3 is the dimensionality, and Vm = cyh? is the “volume” of the marker, where the 
constant cy is straightforward to calculate. 

The above indicates that //^(O) goes from ~ /3/ (ph) for small /? to ~ 1/ {ph) for large (3. At 
intermediate viscous CFL numbers /3, we can set 

//3 (0) = and 5^(0) = 0, 


(18) 
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where C (/3 <C 1) ~ 2/3/(3cy) is linear for small (3 and then becomes 0(1) for large f3. We 
will obtain the actual form of C (/?) from empirical fitting. 


Secondly, for r ^ /i, we know the asymptotic decay of the hydrodynamic interactions from (15): 


• For the steady Stokes regime (/3 —>■ oo), we have the Oseen tensor given by 


focir '> h) ^ goo{r :$> h) ^ (Svrr/r) ^ 


(19) 


• For the inviscid case (/? = 0), we get the electric field of an electric dipole, 

/o(r»/i)«-^^and5o(r-»/i)«^^^, ( 20 ) 

which is also the asymptotic decay for /3 > 0 for r » hy/~p. 

We obtain the actual form of the functions fjsir) and gp(r) empirically by fitting numerical data 
for the parallel and perpendicular mobilities 

4j = PS f0 (nj) + gi3 (r^), 

~ U i^ij) > 

where rjj ■ Vij = 0. To do so, we placed a large number of markers in a cube of length l/S inside 
a periodic domain of length 1. For each marker z, we applied a unit force A* with random direction 
while leaving Aj = 0 for j ^ i, solved ([^, and then interpolated the fluid velocity v at the position 
of each of the markers. The resulting parallel and perpendicular relative velocity for each of the 
N{N — l)/2 pairs of particles allows us to estimate //? (r^) and gp (rij). By making the number 
of markers N sufficiently large, we sample the mobility over essentially all relative positions of the 
pair of markers. For the self-mobility Mu {ru = 0), we take 5 / 3 ( 0 ) = 0 and compute fp{0) from 
the numerical data. 

If the spatial discretization were perfectly translationally and rotationally invariant and the 
domain were infinite, all of the numerical data points for //3 (r) and 5/3 (r) would lie on a smooth 
curve and would not depend on the actual position of the pair of markers relative to the underlying 
grid. In reality, it is not possible to achieve perfect translational invariance with a kernel of finite 
support [ 6 ], and so we expect some (hopefully small) scatter of the points around a smooth fit. 
Normalized numerical data for f/sir) and 5/3 (r) are shown in Fig. and we indeed see that the data 
can be fit well by smooth functions over the whole range of distances. To maximize the quality of the 
fit, we perform separate fits for /3 —)■ 00 (steady Stokes flow) and finite (3. We also make an effort to 
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Figure 2: Normalized mobility functions f{x) (left) and g{x) (right) defined similarly to (All as a function 
of marker-marker distance x = r/h, in three dimensions for the six-point kernel of Bao et al. |48) . over a 
range of viscous CFL numbers (different colors, see legend). Numerical data is shown with symbols and 


obtained using a 256^ periodic fluid grid, while dashed lines show our empirical fit of the form (A2| for 
steady Stokes {jS —)■ oo) and (A6) for finite /3. For steady Stokes flow, the numerical data is in reasonable 
agreement with the RPY tensor (dashed red line). 


make the fits change smoothly as (5 grows towards infinity, as we explain in more detail in Appendix 
Code to evaluate the empirical fits described in Appendices and is publicly available to 
others for a number of kernels constructed by Peskin and coworkers (three-, four-, and six-point) 
in both two and three dimensions at http://cims.nyu.edU/~donev/src/MobilityFunctions.c, 


V. LINEAR SOLVER 

To solve the constrained Stokes problem ([^, we use the preconditioned flexible GMRES (FGM- 
RES) method, which is a Krylov solver. We will refer to this as the “outer” Krylov solver, as it must 
be distinguished from “inner” Krylov solvers used in the preconditioner. Because we use Krylov 
solvers in our preconditioner and because Krylov solvers generally cannot be expressed as linear 
operators, it is crucial to use a flexible Krylov method such as EGMRES for the outer solver. The 
overall method is implemented in the open-source immersed-boundary adaptive mesh refinement 
(IBAMR) software infrastructure |42] : in this work we focus on uniform grids and do not use the 
AMR capabilities of IB AMR (but see O [TT] ). IB AMR uses Krylov solvers that are provided by 
the PETSc library [61] . 
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A. Preconditioner for the constrained Stokes system 

In the preconditioner used by the onter Krylov solver, we want to approximately solve the nested 
saddle-point linear system 


A g -s 


V 


9 

o 

o 

1 


vr 

= 

h 

-J" 0 0 


A 


W 


where we recall that A, = {p/ /S.t)I — ph~‘^Cy. Let us set a = 1 if ^ has a null-space, (e.g., for 
a fully periodic domain for steady Stokes flow) and we set a = 0 if ^ is invertible. When a = 1, 
let ns define the restricted inverse AT^ to only act on vectors of mean value zero, and to return a 
vector of mean zero. 

Applying our Schur complement based preconditioner for solving consists of the following 
steps: 

1. Solve the (unconstrained) fluid sub-problem. 


A 

G 


V 


9 

-T> 

0 


TT 


h 


To control the accuracy of the solution one can either nse a relative tolerance based stopping 
criterion or fix the number of iterations Ng in the inner solver. 

2 . Calculate the slip velocity on the set of markers, AV = — {J'v -|- W). 

3. Approximately solve the Schnr complement system, 

MK = AF, (21) 


4. Optionally, re-solve the corrected fluid sub-problem. 
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-|-«SA — avol ^l^A 

-T> 
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TT 
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where the mobility approximation A4 is constructed as described in Section 


IV 


All linear solvers used in the preconditioner can be approximate, and this is in fact the key to 
the efficiency of the overall solver approach. Notably, the inner Krylov solvers used to solve the 
unconstrained Stokes snb-problems in steps 1 and 4 above can be done by using a small number Ng 
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of iterations using a method briefly described in the next section. If the fluid sub-problem is ap¬ 
proximately solved in both steps 1 and 4, which we term the full Schur complement preconditioner, 


each application of the preconditioner requires 2Ns applications of the Stokes preconditioner (22). 


It is also possible to omit step 4 above to obtain a block lower triangular Schur preconditioner 
|62j . which requires only Ng applications of the unconstrained Stokes preconditioner ( |22[ ). We will 
numerically compare these two preconditioners and study the effect of Ng on the convergence of 
the FGMRES outer solver in Section IVII Al 


B. Unconstrained Fluid Solver 


A key component we rely on is an approximate solver for the unconstrained Stokes sub-problem, 
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TT 
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for which a number of techniques have been developed in the finite-element context [HI]- To 
solve this system, we use GMRES with a preconditioner based on the projection method, 
as proposed by Griffith [32] and improved to some extent by Gai et al. [33]. Specifically, the 
preconditioner for the Stokes system that we use in this work is 

\ / X 0 \ / 0 \ 

( 22 ) 



-1 


--1 


where Cp = {T^G) is the dimensionless pressure (scalar) Laplacian, and A. and Cp denote 

approximate solvers obtained by a single V-cycle of a geometric multigrid solver for the vector 
Helmholtz and scalar Poisson problems, respectively. In the time-dependent case, the approximate 
Schur complement for the unconstrained Stokes sub-problem is 

and for steady Stokes flow, IB = pX. Eurther discussion of the relation of these preconditioners 
to the those described in the book [62| can be found in |32j . 

Observe that one application of is relatively inexpensive and involves only a few scalar 
multigrid V-cycles. Indeed, solving the Stokes system using GMRES with this preconditioner is 
only a few times more expensive than solving a scalar Poisson problem, even in the steady Stokes 
regime jSSj. Note that it is possible to omit the upper right off-diagonal block in the first matrix on 


the right hand side of (22) to obtain a block lower triangular preconditioner that is also effective. 
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and may in fact be preferred at zero Reynolds number since it allows one to skip a sweep of 
the pressure multigrid solver [33]. We empirically find that including the Poisson solve (velocity 
projection) improves the overall performance of the outer solver. 


C. Mobility Solver 


From a computational perspective, one of the most challenging steps in our preconditioner is 


solving the mobility sub-problem (21). Since this is done inside a preconditioner, and because A4 


is itself an approximation of the true mobility matrix A4, it is not necessary to solve (21) exactly. 


In the majority of the examples presented herein, we solve (21) using direct solvers provided by 


LAPACK. This is feasible on present hardware for up to around 10® markers and allows us to focus 
on the design of the approximation A4 and to study the accuracy of the overall method. 

Let us denote with s the smallest marker-marker spacing. For well-spaced markers, s/h ^ 2, 
our approximate mobility A4 is typically SPD even for large numbers of markers, and in these 


cases, we can use the Cholesky factorization to solve (21). In some cases, however, there may 


be a few small or even negative eigenvalues of A4 that have to be handled with care. We have 
found that the most robust (albeit expensive) alternative is to perform an SVD of A4, and to 
use a pseudoinverse of A4 (keeping only eigenvalues larger than some tolerance esvD > 0) to 


solve (21). This effectively filters out the spuriously small or negative eigenvalues. Note that 


the factorization of A4 needs to be performed only once per constrained Stokes solve since the 
body is kept fixed during a time step. In cases where there is a single body, the factorization 
needs to be performed only once per simulation and can be reused; if the body is translating or 


rotating, one ought to perform appropriate rotations of the right hand side and solution of (21). 


In some cases of practical interest where the number of markers is not too large, it is possible 
to precompute the true mobility A4o with periodic boundary conditions (for a sufficienly large 
domain) and to store its factorization. Even if the structure moves relative to the underlying grid, 
such a precomputed (reference) mobility A4o is typically a much better approximation to the true 
mobility than our empirical approximation AL, and can effectively be used in the preconditioner. 
Determining effective approaches to solving the mobility sub-problem in the presence of mnltiple 
moving rigid bodies remains future work, as discussed further in the Conclusions. 
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VI. CONDITIONING OF THE MOBILITY MATRIX 

The conditioning of the constrained Stokes problem Q is directly related to the conditioning 
of the Schur complement mobility matrix A4 = which is intimately connected to the 

relation between the fluid solver grid spacing h and the smallest inter-marker spacing s. Firstly, 
it is obvious that if two markers i and j are very close to each other, then the fluid solver cannot 
really distinguish between Aj and Kj and will instead effectively see only their sum. We also 
know that using too many markers for a fixed fluid grid will ultimately lead to a rank-deficient 
A4, because it is not possible to constrain a finite-dimensional discrete fluid velocity at too many 
points. This physical intuition tells us that the condition number of A4 should increase as the 
marker spacing becomes small compared to the grid spacing. This well-known intuition, however, 
does not tell us how closely the markers can or must be placed in practice. Standard wisdom for 
the immersed boundary method, which is based on the behavior of models of elastic bodies, is to 
make the marker spacing on the order of half a grid spacing. As we show, this leads to extremely 
ill-conditioned mobility matrices for rigid bodies. We note that the specific results depend on the 
dimensionality, the details of the fluid solver, and the specific kernel used; however, the qualitative 
features we report appear to be rather general. 

To determine the condition number of the mobility matrix, we consider “open” and “filled” 
sphere models. We discretize the surface of a sphere as a shell of markers constructed by a recursive 
procedure suggested to us by Charles Peskin (private communication). We start with 12 markers 
placed at the vertices of an icosahedron, which gives a uniform triangulation of a sphere by 20 
triangular faces. Then, we place a new marker at the center of each edge and recursively subdivide 
each triangle into four smaller triangles, projecting the vertices back to the surface of the sphere 
along the way. Each subdivision approximately quadruples the number of vertices, with the fc-th 
subdivision producing a model with 10 • 4^“^ -|- 2 markers. To create filled sphere models, we 
place additional markers at the vertices of a tetrahedral grid filling the sphere that is constructed 
using the TetGen library, starting from the surface triangulation described above. The constructed 
tetrahedral grids are close to uniform, but it is not possible to control the precise marker distances 
in the resulting irregular grid of markers. We use models with approximately equal edges (distances 
between nearest-neighbor markers) of length ~ s, which we take as a measure of the typical marker 
spacing. We numerically computed the mobility matrix A4 for an isolated spherical shell in a large 
periodic domain for various numbers of markers N. Here we keep the ratio s/h fixed and keep 
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the marker spacing fixed at s ~ 1; one can alternatively keep the radius of the sphere fixed In 
Fig. we show the spectrum of for varying levels of resolution for three different spacings of 
the markers, sjh = 1, sjh = 3/2, and sjh = 2. Similar spectra, but with somewhat improved 
condition number (i.e., fewer smaller eigenvalues), are seen for nonzero Reynolds numbers (finite 

/ 3 ). 

The results in Fig. strongly suggest that as the number of markers increases, the low-lying 
(small eigenvalue) spectrum of the mobility matrix approaches a limiting shape. Therefore, the 
nontrivial eigenvalues remain bounded away from zero even as the resolution is increased, which 
implies that for s//i > 1 the system ([^ is uniformly solvable or “stable” under grid rehnement. 
Note that in the case of a sphere, there is a trivial zero eigenvalue in the continuum limit, which 
corresponds to uniform compression of the sphere; this is reflected in the existence of one eigenvalue 
much smaller than the rest in the discrete models. Ignoring the trivial eigenvalue, the condition 
number of M. is 0{N) for this example because the largest eigenvalue in this case increases like the 
number of markers N, in agreement with the fact that the Stokes drag on a sphere scales linearly 
with its radius. This is as close to optimal as possible, because for the continuum equations for 
Stokes flow around a sphere, the eigenvalues corresponding to spherical harmonic modes scale like 
the index of the spherical harmonic. However, what we are concerned here is not so much how the 
condition number scales with N, but with the size of the prefactor, which is determined by the 
smallest nontrivial eigenvalues of Ad. 

Figure clearly shows that the number of very small eigenvalues increases as we bring the 
markers closer to each other, as expected. The increase in the conditioning number is quite rapid, 
and the condition number becomes 0(10® —10^) for marker spacings of about one per fluid grid cell. 
For the conventional choice s ~ h/2, the mobility matrix is so poorly conditioned that we cannot 
solve the constrained Stokes problem in double-precision floating point arithmetic. Of course, if the 
markers are too far apart then fluid will leak through the wall of the structure. We have performed 
a number of heuristic studies of leak through flat and curved rigid walls and concluded that s/h k, 2 
yields both small leak and a good conditioning of the mobillity, at least for the six-point kernel 
used here [IH]. Therefore, unless indicated otherwise, in the remainder of this work, we keep the 
markers about two grid cells apart in both two and three dimensions. It is important to emphasize 
that this is just a heuristic recommendation and not a precise estimate. We remark that Taira and 

^ The scaling used here, keeping s = 1 fixed, is more natural for examining the small eigenvalues of Al, which are 
dominated by discretization effects, as opposed to the large eigenvalues, which correspond to physical modes of 
the Stokes problem posed on a sphere and are insensitive to the discretization details. 
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Figure 3: Eigenvalue spectrum of the mobility matrix for steady Stokes flow around a spherical shell covered 
with different numbers of markers (42, 162, 642, or 2562, see legend) embedded in a periodic domain. Solid 
lines are for marker spacing of s « 2h, dashed-dotted lines for spacing s « 1.5 h, and dashed lines for spacing 
of s « Ih. The marker spacing is s « 1 in all cases; for s « 2h, the fluid grid size is 128^ for 2562 markers 
and 64^ for smaller number of markers, and scaled accordingly for other spacings. For comparison, we show 
the spectrum of M.rpy for the most resolved model {N = 2562 markers) at s/h « 2. Also shown is the 
spectrum of the empirical (fit) approximation to the mobility A4 for the two larger spacings; for s ~ h our 
empirical approximation is very poor and includes many spurious negative eigenvalues (not shown). 

Colonius, who solve a different Schur complement “modified Poisson equation”, recommend s//i ~ 1 
to “achieve a reasonable condition number and to prevent penetration of streamlines caused by a 
lack of Lagrangian points.” 

It is important to observe that putting the markers further than the traditional wisdom will 
increase the “leak” between the markers. For rigid structures, the exact positioning of the markers 
can be controlled since they do not move relative to one another as they do for an elastic bodies; 
this freedom can be used to reduce penetration of the flow into the body by a careful construction 
of the marker grid. In the Conclusions, we discuss alternatives to the traditional marker-based IB 
method [l6] that can be used to control the conditioning number of the Schur complement and 
allow for more tightly-spaced markers. 

It is worthwhile to examine the underlying cause of the ill-conditioning as the markers are 
brought close together. One source of ill-conditioning comes from the discrete (finite-dimensional) 
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nature of the fluid solver, which necessarily limits the rank of the mobility matrix. But another 
contributor to the worsening of the conditioning is the regularization of the delta function. Observe 
that for a true delta function (a —?• 0) in Stokes flow, the pairwise mobility is the length-scale-free 
Oseen tensor ~ and the shape of the spectrum of the mobility matrix has to be independent 
of the spacing among the markers. In the standard immersed boundary method, a ~ /i, so the 
fluid grid scale h and the regularization scale a are difficult to distinguish. 

To try to separate h from a, we can take a continuum model of the fluid, but keep the discrete 
marker representation of the body; see 0. In this case the pairwise mobility would be given by 
(13), which leads to the RPY tensor (|16|) for a kernel that is a surface delta function over a sphere 


of radius a (see (4.1) in [59]). In Fig. [^we compare the spectra of the discrete mobility A4 with 
those of the analytical mobility approximation Adupy constructed by using (16) for the pairwise 
mobility. We observe that the two are very similar for s ~ 2h, however, for smaller spacings 
•^RPY does not have very small eigenvalues and is much better conditioned than A4 (data not 
shown). In Fig. [^we also show the spectrum of our approximate mobility Ad constructed using 
the empirical fits described in Section IV The resulting spectra show a worsening conditioning for 
spacing s ^ 1.5h consistent with the spectrum of Ad. These observations suggest that both the 
regularization of the kernel and the discretization artifacts contribute to the ill-conditioning, and 
suggest that it is worthwhile to explore alternative discrete delta function kernels in the context of 
rigid-body IB methods. 

We also note that we see a severe worsening of the conditioning of A4, independent of j3, when 
we switch from a spherical shell to a filled sphere model. Some of this may be due to the fact that 
the tetrahedral volume mesh used to construct the marker mesh is not as uniform as the surface 
triangular mesh. We suspect, however, that this ill-conditioning is primarily physical rather than 
numerical, and comes from the fact that the present marker model cannot properly distinguish 
between surface tractions and body (volume) stresses. Therefore, A remains physically ill-defined 
even if one gets rid of all discretization artifacts. 

Lastly, it is important to emphasize that in the presence of ill-conditioning, what matters in 
practice are not only the smallest eigenvalues but also their associated eigenvectors. Specifically, 
we expect to see signatures of these eigenvectors (modes) in A, since they will appear with large 
coefficients in the solution of ([^ if the right hand side has a nonzero projection onto the cor¬ 
responding mode. As expected, the small-eigenvalue eigenvectors of the mobility correspond to 
high-frequency (in the spatial sense) modes for the forces A. Therefore, if the markers are too 
closely spaced the solutions for the forces A will develop unphysical high-frequency oscillations or 
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jitter, even for smooth flows, especially in time-dependent flows, as observed in practice [63]. We 
have observed that for smooth flows (i.e., smooth right hand-side of Q), the improved translational 
invariance of the 6-point kernel reduces the magnitude of this jitter compared to the traditional 
Peskin four-point kernel. 


VII. NUMERICAL TESTS 

In this section we apply our rigid-body IB method to a number of benchmark problems. We first 
present tests of the preconditioned FGMRES solver, and then demonstrate the advantage of our 
method over splitting-based direct forcing methods. We further consider a simple test problem at 
zero Reynolds number, involving the flow around a fixed sphere, and study the accuracy of both the 
fluid (Eulerian) variables v and vr, as well as of the body (Lagrangian) surface tractions represented 
by A, as a function of the grid resolution. We finally study flows around arrays of cylinders in two 
dimensions and spheres in three dimensions over a range of Reynolds numbers, and compare our 
results to those obtained by Ladd using the Lattice-Boltzmann method isaissiiMj. 


A. Empirical convergence of GMRES 


Here we consider the model problem of flow past a sphere in a cubic domain that is either periodic 
or with no-slip boundaries. Except for the largest resolutions studied here, the number of markers 


is relatively small, and dense linear algebra can be used to solve the mobility subproblem (21) 


robustly and efficiently, so that the cost of the solver is dominated by the fluid solver. We therefore 


use the number of total applications of the Stokes preconditioner (22) as a proxy for the CPU 


effort, instead of relying on elapsed time, which is both hardware and software dependent. A key 
parameter in our preconditioner is the number of iterations Ng used in the iterative unconstrained 
Stokes solver. We recall that in the full preconditioner, there are two unconstrained inexact Stokes 
solves per iteration, giving a total of 2Ns applications of per outer EGMRES iteration. If the 
lower triangular preconditioner is used, then the second inexact Stokes solve is omitted, and we 
perform only Ng applications of per outer EGMRES iteration. 

In the first set of experiments, we use the full preconditioner and periodic boundary conditions. 
We represent the sphere by a spherical shell of markers that is either empty (162 markers) or is 
filled with additional markers in the interior (239 markers). The top panels of Eig. show the 
relative EGMRES residual as a function of the total number of applications of for several 
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Figure 4: FGMRES convergence for the constrained problem ([^ for different numbers of iterations Ng 
in the unconstrained Stokes solver used in the preconditioner. The specific problem is a rigid sphere of 
hydrodynamic radius R moving through a stationary domain of length L « 4.35 i?, with the marker spacing 
fixed at s/h 2 and the GMRES restart frequency set to 100 iterations. (Top left panel) Steady Stokes 
{Re = 0) flow for an empty 162-marker shell and a filled 239-marker sphere moving through a periodic 
domain of 32^ fluid grid cells. (Top right panel) Same as top left but for Re « 10. (Bottom left panel) 
As top left panel but now in domain with no slip boundary conditions applied on all sides of the domain. 
(Bottom right panel) A spherical shell moving in a non-periodic domain (as in bottom left panel) for different 
resolutions of the shell (162, 642, 2562, and 10242 markers, respectively) and the fluid solver grid (32^, 64^, 
128^, and 256^ grid cells, respectively), fixing Ng = 4, for both the full Schur complement preconditioner 
and the lower triangular approximate Schur complement preconditioner. 


different choices of Ng, for both steady Stokes flow (left panel) and a flow at Reynolds number 
Re = 10 (right panel). We see that for spherical shells with well-conditioned A4 and A4, the 
exact value of Ng does not have a large effect on solver performance. However, making Ng very 
large leads to wasted computational effort by “over-solving” the Stokes system. This degrades the 
overall performance, especially for tight solver tolerance. For the ill-conditioned case of a filled 
sphere model in steady Stokes flow, the exact value of Ng strongly affects the performance, and 
the optimal value is empirically determined to be Ng = 2. As expected, the linear system (§ is 
substantially easier to solve at higher Reynolds numbers, especially for the fllled-sphere models. 

In the bottom left panel of Fig. |^we show the FGMRES convergence for a non-periodic system. 
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In this case, we know that the Stokes preconditioner itself does not perform as well as in the 
periodic case [32l [33] , and we expect slower overall convergence. In this case, we see that Ng = 2 
and = 4 are good choices. Investigations (data not shown) show that = 4 is more robust 
for problems with a larger number of markers. Also, note that increasing Ng decreases the total 
number of FGMRES iterations for a fixed number of applications of 'Pg^, and therefore reduces 


the overall memory usage and the number of times the mobility subproblem (21) needs to be 


solved; however, note that each of these solves is just a backward/forward substitution if a direct 
factorization of A4 has been precomputed. 

The bottom right panel of Fig. shows the FGMRES convergence for a non-periodic system 
as the resolution of the grid and the spherical shell is refined in unison, keeping Ng = 4. The 
results in Fig. demonstrate that our linear solver is able to cope with the increased number 
of degrees of freedom under refinement relatively robustly, although a slow increase of the total 
number of FGMRES iterations is observed. Comparing the full preconditioner with the lower 
triangular preconditioner, we see that the latter is computationally more efficient overall; this is 
in agreement with experience for the unconstrained Stokes system [33]. In some sense, what this 
shows is that it is best to let the FGMRES solver correct the initial unconstrained solution for 
the velocity and pressure in the next FGMRES iteration, rather than to re-solve the fluid problem 
in the preconditioner itself. However, if very tight solver tolerance is required, we find that it is 
necessary to perform some corrections of the velocity and pressure inside the preconditioner. In 
principle, the second unconstrained Stokes solve in the preconditioner can use a different number 
of iterations Ng from the first, but we do not explore this option further here. Also note that if 
Ad ~ Ad (for example, if it was computed numerically rather than approximated), then the full 
Schur complement preconditioner will converge in one or two iterations and there is no advantage 
to using the lower triangular preconditioner. 


B. Flow through a nozzle 


In this section we demonstrate the strengths of our method on a test problem involving steady- 
state flow through a nozzle in two dimensions. We compare the steady state flow through the 
nozzle obtained using our rigid-body IB method to the flow obtained by using a splitting-based 
direct forcing approach mm- Specifically, we contrast our monolithic fluid-solid solver to a split 
solver based on performing the following operations at time step n:Solve the fluid sub-problem as 
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1. Calculate the slip velocity on the set of markers, AV = — giving 

the fluid-solid force estimate = {p/At) AV. 

2. Correct the fluid velocity to approximately enforce the no-slip condition, 

^n+l ^ AV. 

Note that in the original method of [2] in the last step the fluid velocity is projected onto the space of 
divergence-free vector fields by re-solving the fluid problem with the approximation {p/ At) I 
(i.e., ignoring viscosity). We simplify this step here because we have found the projection to make 
a small difference in practice for steady state flows, since the same projection is carried out in the 
subsequent time step. 

We discretize a nozzle constriction in a slit channel using IB marker points about 2 grid spacings 
apart. The geometry of the problem is illustrated in the top panel of Fig. parameters are p = 1, 
grid spacing Ax = 0.5, nozzle length I = 55.5, nozzle opening width d ~ 2.9, and rj variable (other 
parameters are given in the figure caption). No slip boundary conditions are specified on the top 
and bottom channel walls, and on the side walls the tangential velocity is set to zero and the normal 
stress is specified to give a desired pressure jump across the channel of Avr = 2. The domain is 
discretized using a grid of 256 x 128 cells and the problem evolved for some time until the flow 
becomes essentially steady. The Reynolds number is estimated based on the maximum velocity 
through the nozzle opening and the width of the opening. 

In the bottom four panels in Fig. [^we compare the flow computed using our method (left panels) 
to that obtained using the splitting-based direct forcing algorithm summarized above (right panels). 
Our method is considerably slower (by at least an order of magnitude) for this specific example 
because the GMRES convergence is slow for this challenging choice of boundary conditions at 
small Reynolds numbers in two dimensional (recall that steady Stokes flow in two dimensions has 
a diverging Green’s function). To make the comparison fairer, we use a considerably smaller time 
step size for the splitting method, so that we approximately matched the total execution time 
between the two methods. Note that for steady-state problems like this one with fixed boundaries, 
it is much more efficient to precompute the actual mobility matrix (Schur complement) once at the 
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Figure 5: Two-dimensional flow through a nozzle (top panel) in a slit channel computed by our rigid IB 
method (left panels) and a simplified version of the splitting-based method of Bhalla et al. [2] (right panels). 
(Top panel) The geometry of the channel along with the (approximately) steady state flow at Rewlfl, as 
obtained using our method. The color plot shows the pressure and the velocity is shown as a vector held. 
(Middle panels) Flow at Re « 19, computed at time T = 10^. For our method (left) we use a time step size 
of At = 5 • 10“^ (corresponding to advective CFL number of Umax^t/Ax « 0.13), while for the splitting 
method (right) we use At = 10“^. The streamlines are traced from the entrance to the channel for a time 
of Tg = 7 • 10^ and shown as black lines. (Bottom panels) Same as the middle row but now for Re « 0.2, 
final time T = 10 and streamlines followed up to Tg = 4 • 10"^, with At = 0.125 for our method (left), and 
At = 10“^ for the splitting method (right). 
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beginning, instead of approximating it with our empirical fits. However, for a more fair and general 
comparison we instead use our preconditioner to solve the constrained fluid problem in each time 
step anew to a tight GMRES tolerance of 10“®. For this test we use Ng = 2 iterations in the fluid 
solves inside our preconditioner. 

The visual results in the right panels of Fig. clearly show that the splitting errors in the 
enforcement of the no-slip boundary condition lead to a notable “leak” through the boundary, 
especially at small Reynolds numbers. To quantify the amount of leak we compute the ratio of the 
total flow through the opening of the nozzle to the total inflow; if there is no leak this ratio should 
be unity. Indeed, this ratio is larger than 0.99 for our method at all Reynolds numbers, as seen in 
the lack of penetration of the flow inside the body in the left panels in Fig. For Re ~ 19, we 
find that even after reducing the time step by a factor of 50, the splitting method gives a ratio of 
0.935 (i.e., 6.5% leak), which can be seen as a mild penetration of the flow into the body in the 
middle right panel in Fig. For Re ~ 0.2, we find that we need to reduce At by a factor of 1250 
to get a flow ratio of 0.94 for the splitting method; for a time step reduced by a factor of 125 there 
is a strong penetration of the flow through the nozzle, as seen in the bottom right panel of Fig. 

C. Stokes flow between two concentric shells 

Steady Stokes flow around a fixed sphere of radius Ri in an unbounded domain (with fluid at 
rest at infinity) is one of the fundamental problems in fluid mechanics, and analytical solutions 
are well known. Our numerical method uses a regular grid for the fluid solver, however, and thus 
requires a finite truncation of the domain. Inspired by the work of Balboa-Usabiaga et al. [65], 
we enclose the sphere inside a rigid spherical shell of radius R 2 = 4Ri. This naturally provides a 
truncation of the domain because the flow exterior to the outer shell does not affect the flow inside 
the shell. Analytical solutions remain simple to compute and are given in Appendix [C| 

We discretize the inner sphere using a spherical shell of markers, since for steady Stokes flow 
imposing a rigid body motion on the surface of the inner sphere guarantees a stress-free rigid body 
motion for the fluid filling the inner sphere |2T]. We use the same recursive triangulation of the 
sphere, described in Section |V| to construct the marker grid for both the inner and outer shells. 
The ratio of the number of markers on the outer and inner spheres is approximately 16 (i.e., there 
are two levels of recursive refinement between the inner and outer shells), consistent with keeping 
the marker spacing similar for the two shells and a fixed ratio R 2 /R 1 = 4. The fluid grid size 
is set to keep the markers about two grid cells apart, s k. 2h. The rigid-body velocity is set to 
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Figure 6: Marker configuration for computing Stokes flow between two concentric spherical shells. Markers 
are shown as spheres with size on the order of their effective hydrodynamic radius. The inner shell of markers 
is shown in red, and the outer shell of markers is shown in gray. (Left) Intermediate resolution, inner shell 
of 162 markers and outside shell of 2562 markers. (Right) Highest resolution studied here, inner shell of 642 
markers and outside shell of 10242 markers. 

V = (1,0,0) for all markers on the outer shell, and to V = 0 on all markers on the inner shell. The 
outer sphere is placed in a cubic box of length I = 4.15i?2 with specihed velocity v = (1,0,0) on all 
of the boundaries; this choice ensures that the flow outside of the outer shell is nearly uniform and 
equal to u = (1,0,0). In the continuum setting, this exterior flow does not affect the flow of interest 
(which is the flow in-between the two shells), but this is not the case for the IB discretization since 
the regularized delta function extends a few grid cells on both sides of the spherical shell. 

A spherical shell of geometric radius Rg covered by markers acts hydrodynamically as a rigid 
sphere of effective hydrodynamic radius Rh ^ Rg + a [66], where a is the hydrodynamic radius 
of a single marker na Ea EH] (we recall that for the six-point kernel used here, a ~ 1.47 h). 
A similar effect appears in the Lattice-Boltzmann simulations of Ladd, with a being related to 
the lattice spacing isiiEa 166] . When comparing to theoretical expressions, we use the effective 
hydrodynamic radii of the spherical shells (computed as explained below) and not the geometric 
radii. Of course, the enhancement of the effective hydrodynamic radius over the geometric one 
is a numerical discretization artifact, and one could choose not to correct the geometric radius. 
However, this comparison makes immersed boundary models of steady Stokes flow appear much 
less accurate than they actually are in practice. For example, one should not treat a line of markers 
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as a zero-thickness object of zero geometric radius; rather, such a line of rigidly-connected markers 
should be considered to model a rigid cylinder with finite thickness proportional to a |12j . 

We can measure the effective hydrodynamic radius Rh of a spherical shell of markers from the 
drag force on a periodic cubic lattice of such objects moving with velocity V. Specifically, we 
place a single shell of N markers in a triply-periodic domain with cubic unit cell of length I, set 
V = (1,0,0) on all markers, solve and measure the total drag force as F = 
periodic correction to the Stokes drag formula is well-known 163, 


F 


67rRh 


(23) 


T]V 1 - 2.8373{Rh/l) + A.19{Rh/l)^ - 27A{Rh/l)^ + h.o.t. ’ 
and allows us to obtain a very accurate estimate of Rh from the drag for I ^ Rh- The results are 
given in the left half of table [l] in the form of the dimensionless ratio Rh/Rg', we see that as the 
resolution is increased Rh —?■ Rg with an approximately linear rate of convergence, as expected. 
Since this computation refers to flow outside of the shell of markers, we can call the computed Rh 
the effective outer hydrodynamic radius and use it to set the value of Ri in the theory. We use 
a similar procedure to measure an effective inner hydrodynamic radius i?2 for the outer spherical 
shell. Specifically, we obtain R 2 from the drag on the inner sphere based on the theoretical formula 


(Cl), where we use the previously-determined value of Ri for the effective radius of the inner 
sphere. The results are given in the right half of Table |I] and again show that as the grid is refined 
the hydrodynamic radii converge to the geometric ones. 


Table I: Ratio of the effective hydrodynamic and geometric radii of the inner (left half) and outer (right 
half) spherical shells for simulations of steady Stokes flow around a fixed sphere embedded within a moving 
spherical cavity, at different resolutions. 


Resolution 

Number markers 

Rh/Rg 

Number markers 

Rh/Rg 

grid size 

Inner shell 


Outer shell 


30^ 

12 

1.48 

162 

0.93 

60^ 

42 

1.22 

642 

0.96 

120^ 

162 

1.09 

2562 

0.98 

240^ 

642 

1.04 

10242 

0.99 


1. Convergence of fluid flow (pressure and velocity) 


The top panel of Fig. shows a slice through the middle of the nested spherical shells along 
with the fluid velocity v. Recall that the flow inside the inner sphere should vanish, implying that 
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the pressure inside the inner shell should be constant (set to zero here), and the flow outside of the 
outer sphere should be uniform. The bottom right panel of the figure zooms in around the inner 
sphere to reveal that there is some spurious pressure gradient and an associated counter-rotating 
vortex flow generated inside the inner sphere. The bottom right panel shows the error in the 
computed fluid flow (u,7r), that is, the difference between the computed flow and the theoretical 
solution given in Appendix [C) It is clear that the majority of the error is localized in the vicinity 
of the inner shell and in the interior of the inner sphere. Note that these errors would be much 
larger if the theory had used the geometric radii instead of the hydrodynamic radii for the shells. 



Figure 7: Flow field around a fixed sphere inside a moving spherical cavity. The outer shell is discretized using 
2562 markers, while the inner one has 162 markers, shown as black circles. (Top) Velocity field. (Bottom 
Left) Zoom of the velocity (vector field) and pressure (color plot) around the inner shell. (Bottom Right) 
Same as bottom left but now showing the error in the velocity and pressure compared to the theoretical 
expressions. 
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Tables and III show the norms of the error in the computed flow held as a function of 
resolution. Asymptotically hrst-order convergence is observed in the Li and L 2 norms for both 
the velocity and the pressure. In the L^o norm, we expect the velocity to also converge linearly, 
but we do not expect to see convergence in the pressure, since the velocity is continuous across the 
interface but the pressure has a jump; this is consistent with the numerical data. 


Table II: Normalized norms of the error in the computed velocity (An) for steady Stokes flow around a fixed 
sphere embedded within a moving spherical cavity, at different resolutions (see two left-most columns). An 
estimated order of convergence based on successive refinements is indicated in the column to the right of 
the corresponding error norm. 


Markers 

Resolution 

l|An||i/||n||i 

Rate 

IIAnIb/llnIb 

Rate 

||An||oo/||n||oo 

Rate 

162-12 

30^ 

4.08 • 10-2 


6.39 • 10-2 


0.558 


642-42 

603 

1.14-10-2 

1.83 

2.08 • 10-2 

1.62 

0.322 

0.79 

2562-162 

1203 

4.61 • 10-3 

1.30 

8.74 • 10-3 

1.24 

0.160 

1.01 

10242-642 

2403 

2.16 • 10-3 

1.09 

4.26 • 10-3 

1.04 

0.091 

0.82 


Table III: Normalized norms of the error in the pressure (Att) for steady Stokes flow around a fixed sphere 
embedded within a moving spherical cavity, at different resolutions (see two left-most columns). An esti¬ 
mated order of convergence based on successive refinements is indicated in the column to the right of the 
corresponding error norm. 


Markers 

Resolution 

l|A^||i/IM|i 

Rate 

||A7r||2/||7r||2 

Rate 

IIAttI loo/I kl loo 

Rate 

162-12 

CO 

0 

CO 

0.849 


0.788 


1.0 


642-42 

CO 

0 

0.567 

0.58 

0.486 

0.70 

1.0 

0.0 

2562-162 

1203 

0.344 

0.72 

0.275 

0.82 

0.860 

0.22 

10242-642 

2403 

0.196 

0.81 

0.164 

0.75 

0.704 

0.29 


2. Convergence of Lagrangian forces (surface stresses) 

The hrst-order convergence of the pressure and velocity is expected and well-known in the im¬ 
mersed boundary community. The convergence of the tractions (cr ■ n) on the huid-body interface 
is much less well studied, however. This is in part because in penalty-based or splitting methods, 
it is difficult to estimate tractions precisely (e.g., for penalty methods using stiff springs, the spring 
tensions oscillate with time), and in part because a large number of other studies have placed the 
markers too closely to obtain a well-conditioned mobility matrix and thus to obtain accurate forces. 
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Furthermore, there are at least two ways to estimate surface tractions in IB methods, as discussed 
in extensive detail in Ref. [68]. One method is to estimate fluid stress from the fluid flow and 
extrapolate toward the boundary. Another method, which we use here, is to use the computed 
surface forces A to estimate the tractions. 



Figure 8: Convergence of surface stresses to their theoretical values for the three different resolutions. 
Pointwise traction estimates are shown with symbols as a function of the angle 9 relative to the direction of 
the flow, while the theory is shown with a solid black line. {Left columns) Normal component of the traction 
an = r ■ cr ■ n. {Center columns) Tangential component of traction in direction of flow, ag = 6 ■ cr ■ n. {Right 
columns) Tangential component in the direction perpendicular to the flow, a,f, = <p ■ cr ■ n, which should 
vanish by symmetry. ( Top row) Different resolutions (see legend) for a fixed spacing s « 2h. Note that 
for the coarsest resolution of only 12 markers on the inner sphere, the computed tractions have values off 
the scale of this plot and are thus not shown. {Bottom row) The most resolved case of 2562 — 162 markers 
for different spacing between markers, as indicated in the legend. Note that using s ks h leads to severe 
ill-conditioning and the computed tractions show random scatter well beyond the scale of the plot and are 
thus not shown. 

We obtain pointwise estimates of the tractions at the positions of the markers from the relation 
{cr ■ n) {Ri) ~ where AAj is the surface area associated with marker i. We obtain AAj 

from the surface triangulation used to construct the marker of grids by assigning one third of the 
area of each triangle to each of its nodes. In Fig. [^we show the computed normal and tangential 
components of the traction in polar coordinates, with the z symmetry axes along the direction of 
the flow. The theoretical prediction given in Appendix is shown with a black line and is based 
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on the geometric radii. 

The top row of Fig. shows the computed tractions for several resolutions with marker spacing 
s ^ 2h. It is seen that as the grid is refined, the computed tractions appear to converge pointwise to 
the correct values. However, the convergence is very slow, and even for the large resolutions reported 
here, it is evident that the asymptotic convergence regime has not been reached. Consequently, no 
precise statement about the order of convergence can be made from these data. At lower resolution, 
some of the results even show qualitatively wrong behavior. For example, the normal traction 
an = n ■ a ■ n for a resolution of 42 inner and 642 outer markers grows with 6, but the theoretical 
result decreases with 6. We also see scatter in the values among individual markers, indicating 
that the geometrical and topological non-uniformity of the marker grid affect the pointwise values. 

Nonetheless, we remark that low-order moments of the surface tractions are much more accurate 
than the pointwise tractions. For example, the total drag on the inner sphere is much more accurate, 
as seen in Table |IJ Other test problems not reported here indicate that stresslets are also computed 
quite accurately, especially if one accounts for the distinction between geometric and hydrodynamic 
radii. These findings suggest that weak convergence of the tractions is more robust than strong 
convergence. In fact, lower order moments can show reasonable behavior even if the marker spacing 
is small and the pointwise forces are numerically unstable to compute. Somewhat unsurprisingly, 
we find that the pointwise traction estimates are improved as the spacing among the markers is 
increased; see the bottom row in Fig. The improvement is not only due to the reduction of the 
scatter, as expected from the improvement in conditioning number of the mobility matrix, but also 
due to global reduction of the error in the tractions; the observed global reduction may, however, 
be specific to steady Stokes flow. For widely spaced markers, however, the error in the computed 
flow field field will increase because the flow will penetrate the shell boundary. This once again 
demonstrates the delicate balance that is required in choosing the marker spacing for rigid bodies, 
as we discuss further in the Conclusions. 

D. Steady Stokes flow around sphere in a slit channel 

In this section we study a problem at zero Reynolds numbers with nontrivial boundary condi¬ 
tions, namely, steady Stokes flow around a sphere in a slit channel (flow between two parallel walls). 
It is well-known that computing flows in such geometries using Green’s function based methods 
such as boundary-integral methods is highly-nontrivial [SZlEaiTni. Specific methods for spheres in 
a channel have been developed m but these are not general, in particular, flow in a square channel 
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requires a different method, and incorporating the periodicity in some of the dimensions is non¬ 
trivial [701E2]. At the same time, we wish to point out that boundary-integral methods have some 
advantages over our IB method as well. Notably, they are considerably more accurate, and han¬ 
dling domains unbounded in one or more directions is possible by using the appropriately-decaying 
Green’s function. 

Unlike the case of a single no-slip boundary, writing down an analytical solution for slit channels 
is complex and requires numerically-evaluating the coefficients in certain series expansions m- For 
the component of the mobility /r = F/V of a sphere in an infinite slit channel, Faxen has obtained 
exact series expansions for the mobility at the half and quarter channel locations. 






^ 

2 / 67rr]Rh 

H=d\ ^ 1 

4/ GTrrjRh 


Rh Rl Ri Rl 

1 - 1 . 004 ^ + 0 . 418 ^ + 0 . 21 ^ - 0 . 169 ^ + ... 
H 


Ri 


Ri 


Rl 


1 - 0.6526^ -b 0.1475^ - 0.131^ - 0.0644^ -b . 
H 


(24) 


where Rh is the (hydrodynamic) radius of the sphere, H is the distance from the center of the 
particle to the nearest wall, and d is the distance between the walls. 

To simulate a spherical particle in a slit channel we place a single spherical shell with different 
number of IB markers in a domain of size L x L x d, at either a quarter or half distance from the 
channel wall. No slip walls are placed at z = 0 and z = L, and periodic boundary conditions are 
applied in the x and y directions. For each L, we compute an effective hydrodynamic radius Rl 
by assuming (24) holds with Rh replaced by Rl- We know that as L —)■ oo we have Rl Rh, 


however, we are not aware of theoretical results for the dependence Rl{L) at finite L. In Fig. we 
plot RL/Rh — 1 versus Rh/L for d^. 8Rh- Here the effective hydrodynamic radius of the shell Rh 
is estimated by using (23), as shown in Table [l] (see inner radius). We see that we have consistent 
data for Rl{L) among different resolutions, and we obtain consistency in the limit L —>■ oo. This 
indicates that even a low-resolution model with as few as 12 markers offers a reasonably-accurate 
model of a sphere of effective radius Rh, independent of the boundary conditions. 


E. Steady Stokes flow around cylinders 


Here we study the drag force on a periodic square array of cylinders (i.e., disks in two dimensions) 
with lattice spacing 1. The corresponding study in three dimensions is presented in Section [VH G[ 


The equivalent of (23) in two dimensions for dilute systems is mm 
F Att 

yV ~ - In(V^) - 0.738 + cj) - 0.88702 + 2.038(/)3 -b O ((/>4) ’ 


( 25 ) 
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Figure 9: Effective hydrodynamic radius Rl{L) of a sphere of hydrodynamic radius Rh, translating parallel 
to the walls of a slit channel of dimensions L x L x d. Red symbols are for the sphere at the midplane of 
the channel, H/d = 0.5, and blue symbols are for H/d = 0.25; these two give different dependence Rl{L) 
as expected. The channel width is taken d « ^Rh and different numbers of markers are used for the sphere 
(see legend), and the grid spacing is set to give a marker spacing s as close as possible to a/s « 0.5. Note 
that similar to the example of flow between concentric spheres the correct value of the drag is determined 
by the larger hydrodynamic and not by the geometric radius of the shell. 


where (j) = is the packing fraction of the disks and is the hydrodynamic radius of 


the cylinder, which is defined from (25). Observe that in two dimensions, there is no limit as 
(/> —7- 0, in agreement with Stokes’s paradox for flow around a single cylinder; one must account 
for inertial effects for very small volume fractions in order to obtain physically-relevant results. 


Table IV reports Rh for several different marker models of a cylinder, as estimated by computing 


the drag for a range of packing fractions and extrapolating to (/> <C 1 using (25). As expected, 
the more resolved the cylinder, the closer Rh is to Rg. Filling the cylinder with markers both 
substantially enlarges the effective hydrodynamic radius and also degrades the conditioning of the 
mobility matrix, and is therefore not advised at zero Reynolds number. 

Another interesting limit for which there are theoretical results is the dense limit, in which the 
disks/cylinders almost touch, so that there is a lubrication flow between them. In this limit 

F Ott _5 

vV 2i 


( 26 ) 
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Table IV: Hydrodynamic radii of several discretizations of a cylinder with different numbers of markers on 
the surface and the interior of the body, keeping s//i « 2. Two models have markers only on the surface of 


the cylinder (see the right panel of Fig. 12 for a 39-marker shell). The rest of the models are constructed 
from a regular polar grid of markers filling the interior of the cylinder (see, for example, left panel of Fig. 


12 for a 121-marker cylinder). 


Number markers 

Surface markers 

Interior markers 

Rh/Rg 

39 shell 

39 

0 

1.04 

121 cylinder 

37 

84 

1.15 

100 shell 

100 

0 

1.02 

834 cylinder 

100 

734 

1.03 


where e = 1 — ^= {I — 2Rh) jl is the relative gap between the particles. Note that because 
the number of hydrodynamic cells must be an integer, we cannot get an arbitrary gap between the 
cylinders for a given cylinder model and fixed sjh (i.e., a fixed Rh/Rg)- Also note that when the 
gap between the cylinders is too small, the kernels from markers on two cylinders start to overlap, 
and the problem becomes ill-conditioned; we have been able to compute reliable results down to a 
relative gap of e > 10“^ for the resolutions studied here. 
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Figure 10: The drag coefficient for a periodic array of cylinders in steady Stokes flow for different resolutions 


(see Table IV). (Left) As function of volume fraction, compared to the results of a highly-accurate boundary- 
integral method. (Right) Zoom in for close-packed arrays with inter-particle gap plotted on a log scale to 
show the asymptotic divergence of the lubrication force. 


Numerical results for the normalized drag over a broad range of volume fractions are shown 


in Fig. 10 and compared to results obtained using an in-house two-dimensional version of the 


spectrally-accurate boundary integral method proposed in Ref. m- We obtain very good agree¬ 
ment, similar to that observed using the Lattice Boltzmann method [52], indicating that even 
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moderately-resolved cylinders are good representations so long as one uses their hydrodynamic 
rather than their geometric radius when computing the effective volume fraction. In particular, in 


the right panel, we obtain excellent agreement with the lubrication result (26), seeing an increase 
in the drag of over six orders of magnitude consistent with theory. Of course, the IB method results 
for the drag do not have a true divergence as e —?• 0 because of the regularization of the singular 
kernel; one must use adaptively-refined non-regularized boundary integral methods to truly resolve 
the divergence. In practice, however, effects not included in the theoretical model, such as surface 
roughness or partial slip, will mollify the unphysical divergence. 


F. Unsteady flow around cylinders 


Next, we examine the ability of our rigid-body IB method to model unsteady two-dimensional 
flow around cylinders (disks). We define Reynolds number by 

pVRh VRh 


Re = 


r? 


where Rh is the hydrodynamic radius of the cylinder measured at Re = 0 (see table IV), and V is 
the velocity of the incident flow. For small Reynolds numbers, the mean drag per unit length F is 
given by [S2] 

R 9 

— = ko + k 2 Re^, 
r/V 

where ko {4>) and k 2 {4>) are constants that depend on the packing fraction cj) (defined using the 
hydrodynamic radius). In the range Re ~ 2 — 5, the drag becomes quadratic in the flow rate 
and for moderate Reynolds numbers, a drag coefficient is defined from the empirical relation 


pvmH 

As the Reynolds number is increased, the flow becomes unsteady and vortex shedding occurs, and 
eventually there is a transition to three-dimensional flow. Here we focus on steady flow at Re < 100. 

A staggered-grid variant of the piecewise-parabolic Godunov method is used for spatial dis¬ 
cretization of the advective terms, as explained in detail by Griffith [32]. In our tests, the time step 
size is determined by fixing the advective Gourant number VAt/h = 0.1; this value is well bellow 
the stability limit and ensures that the discretization errors coming from the (unconstrained) fluid 
solver are small. The Adams-Bashforth method is used to handle advection explicitly. The viscous 
terms are handled implicitly using the backward Euler method rather than the implicit midpoint 
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rule because we are interested in steady states and not transient phenomena. We initialize the 
simulations with the fluid moving at a uniform velocity but allow enough time for a steady-state 
to be reached. 


1. Drag on periodic array of disks 


The permeability of a periodic array of aligned cylinders is a well-studied problem and can be 
computed by placing a single cylinder in a periodic domain. To create flow through the periodic 
system, we follow Ladd et al. [SUES] and apply a constant body force f throughout the domain 
(including in the interior of the body). We solve the constrained time-dependent problem to 
a steady state, keeping the cylinder at rest, and measure the average velocity in the domain, 
V = vol~^ JyVdr. In two spatial dimensions, the dimensionless drag coefficient is defined by 


Ie. 

riVx 


where the force F = vol/ = —l^A is the total force applied to the fluid, which must also equal 
the negative of the total force exerted on the rigid body. 

Theory suggests that the correction to the drag scales as Re^ for small Reyonds numbers due 
to the anti-symmetry of the correction to the flow (relative to steady Stokes) of order Re [S2], so 
that 


k = ko + A:2Re^, (27) 

where the values A:o(</>) and k2{4>) depend on the packing fraction (p. To obtain ko, we move the body 
at a constant velocity and obtain the drag force l^A from the solution of the constrained steady 
Stokes problem 0- Because marker-based models of rigid bodies do not have perfect symmetry, 
the force /g = —vol~^ small nonzero components in the direction perpendicular to the 

flow. To ensure that in the limit Re —?• O"*" we have perfect consistency between the finite Re and 
zero Re computations, we use the force / = (/c//co)/q to drive the flow at finite Re numbers. 
Note that it can take thousands of time steps for the steady state to be established for Re > 1; 
to accelerate convergence, we initialize the computation for a given Re from the steady state for 
the closest smaller Re. Also note that the exact mobility matrix A4 and its factorization can be 
precomputed once at the beginning and used repeatedly for these steady-state calculations. 

Fig. shows the dimensionless excess drag k — ko as a function of Re at packing fraction 

(p = 0.193, which is close to the packing fraction cp = 0.2 studied using the Lattice-Boltzmann 
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Figure 11: The dimensionless excess (over Stokes flow) drag coefficient for a square array of disks with 
packing fraction (j) « 0.193 (129-marker filled cylinder model, fluid grid of 64^ cells). Comparison is made 
to known small-Re dependence of the form fcg + k 2 Re^, with the coefficients ko and k 2 taken from the work 
of Koch and Ladd [52] at 0 = 0.2. 


method in Ref. 


. We see very good agreement of the theoretical formula (27) with our results 


using the values of /cq = 49.2 and k 2 = 0.24, which are in good agreement with the values of 
feo = 51.2 and k 2 = 0.26 given in the caption of Fig. 1 in Ref. 


2. Flow past a periodic column of cylinders 

Here we compute several solutions for flow past a column of cylinders at somewhat larger 
Reynolds numbers, mimicking the setup of Ladd |53) . The domain is a long narrow channel of 
2048 X 128 grid cells ® with grid spacing h = 0.5, keeping the markers at a distance 2h. Periodic 
boundary conditions are used in the direction of the short side of the channel (y). The flow is 
driven by “uniform” inflow and outflow boundary conditions in the long direction (x). Specifically, 
we impose a specified normal velocity V and zero tangential velocity at both ends of the channel 

® For very elongated domains, our multigrid-based preconditioner converges much faster for grid sizes that are powers 
of two. 
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The center of the cylinder is fixed at a quarter channel length from the inlet. The cylinders in 
the periodic column are separated by approximately 10 hydrodynamic radii in the y direction (the 
separation is 9.958Rh for the 121-marker cylinder, and 9.875Rh for the 39-marker shell). 



Figure 12: Steady incompressible flow at Re = 10 past a periodic column of cylinders represented as either 
filled disks of 121 markers (left) or a shell of 39 markers (right). We show the magnitude of the Eulerian 
constraint force <SA (gray color map), the streamlines outside of the wake (solid blue lines), the wake velocity 
field (black arrows), and the Lagrangian constraint forces associated with each marker (color arrows). The 
red arrow marks the stagnation point where Vx = 0, as used to determine the wake length. 


Representative flow fields are shown in Fig. 12 for Re = 10 for a filled cylinder model (left) and 
an empty shell model (right). Note that for the computation of total drag on a fixed cylinder either 
model can be used since the spurious flow seen inside the empty shell does not generate any overall 
acceleration of the fluid inside the body. Also note that the spurious counter-rotating vortex pair 
inside the shell diminishes under refinement, at an approximately linear convergence rate, just as 
for steady Stokes flow. Computed drag coefficients k and wake lengths are shown in Table [V] and 
good agreement is seen with the results of Lattice-Boltzmann and finite difference schemes m- 
The wake length measures the distance from the cylinder center to the stagnation point, which is 
obtained by hnding the largest x coordinate on the contour of zero horizontal velocity, Vx = 0. 


G. Flow past periodic arrays of spheres 

Finally, we study the drag on a cubic arrays of spheres of radius a at zero and finite Reynolds 
numbers, and compare our results to those of Hill et al. m- At small packing (volume) fractions 
(j) and Reynolds numbers, according to Eqs. (1-2) in [6l], F — Fq = 3Re/8-|-h.o.t. if <C Re <C 1, 

® An alternative is to use zero tangential stress on both boundaries, or zero normal and tangential stress on the 
outflow; such stress boundary conditions are supported in the fluid solver in the IBAMR library | 32 |. 
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Table V: Numerical results for steady flow past a periodic column of cylinders at different Reynolds number, 


for two different models of the body (see Fig. 12), either a filled cylinder or an empty shell of markers. For 
comparison we reproduce the results in Table 5 in |53j . which are computed either using either a Lattice- 
Boltzmann (LB) or a finite difference (FD) method. (Left) Mean drag coefficient. (Right) Wake length in 


units of Rh- 


Re 

121 cyl 

39 shell 

LB 

FD 

5 

1.52 

1.40 

1.5 

1.49 

10 

2.55 

2.59 

2.6 

2.65 

20 

4.50 

4.61 

4.7 

4.74 

50 

9.96 

9.91 

10.7 

10.3 


Re 

121 cyl 

39 shell 

LB 

FD 

5 

4.31 

4.35 

4.21 

4.32 

10 

2.96 

2.99 

2.91 

2.98 

20 

2.16 

2.19 

2.17 

2.19 

50 

1.55 

1.58 

1.67 

1.61 


or, more relevant to our study, F — Fq ~ Re^/y/^ if Re <C y/cf) <C 1. For small Re and at larger 
densities, the theoretical arguments in [SD E2] predict that the dimensionless drag is quadratic in 
Re because the linear term vanishes by symmetry, so that 

F 


k = 


ko + A:2Re . 


GTrrjaV 

For larger Re, the dependence is expected to switch to linear in Re. 

Here we focus on close-packed cubic lattices of spheres with packing fraction ((> = 7r/6 ~ 0.5236. 
Note that unlike the case of two spatial dimensions, in three dimensions the flow does not need to 
squeeze in-between the (nearly) touching bodies, so the drag does not diverge even at close packing. 
The value of the steady Stokes drag /cq is tabulated in Table |V^ for several resolutions. Different 
resolutions are examined: an empty shell (see Table of 162 (grid size is 16^) or 642 markers (30^ 
grid), as well as a filled sphere of 56 (42 on the surface, 10^ grid) or 239 (162 on the surface, 16^ 
grid) markers; the actual value of the packing fraction based on the effective hydrodynamic radius 
of the model is indicated in the table. A large difference is seen between the filled and empty 
shell models at this high packing fractions because the spheres are very close to each other and 
discretization artifacts become pronounced. We have also performed simulations at a lower (but 
still high) packing fraction of 4> = 0.44, and there we see much better agreement between the filled 
and empty sphere models; note that at small (/> <C 1 the value of ko must match among resolutions 
since we define the packing fraction from R^, which is itself determined from the value of ko at 


small (f) using (25). 

Numerical results for the dimensionless drag coefficient k near the close-packed density ~ 0.52 


are shown in Fig. 13 Because our discrete models of spheres do not have the same symmetry as 


a perfect sphere, we numerically observe a small 0(Re) correction that can dominate the true 
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Table VI: Dimensionless drag force kg for steady Stokes flow {Re = 0) past a simple-cubic array of spheres 
at volume fraction « 7r/6 (close packing). For the highest-resolution LB simulations in the reported 
value is kg = 42.8. 


Number of markers 

(t> 

fco 

56 filled 

0.5236 

40.08 

239 filled 

0.5238 

40.73 

162 shell 

0.5213 

44.49 

642 shell 

0.5236 

43.29 


correction A: 2 Re^ for Re <C 1; this is especially evident in the right panel of Fig. 13 for coarsely- 
resolved models (e.g., a 56-marker sphere) Empirical fits to literature data for ko, ki, k 2 , and the 


range of Re values over which the various fits are valid are tabulated in Ref. m- Fig. [l^ compares 
our results to these fits, as well as to reference results obtained using the Lattice Boltzmann method 
[53] at close packing. We observe the expected switch from linear to quadratic dependence on Re 
and also a reasonable agreement with the literature data, and the agreement appears to improve 
with increasing resolution. 



Figure 13: Numerical values (symbols) for the drag coefficient of a periodic array of spheres close packed in 
a cubic lattice of volume fraction (f) = 7r/6 « 0.52, for several resolutions (see legend), using a linear (left) 
or log scaling (right). Comparison is made to empirical formulas given in |73j (lines), as well as Lattice 
Boltzmann results for close-packed cubic arrays given in Table 6 in [SS] (crosses). 


In two dimensions, we can more easily make the discrete models symmetric and this is why Fig. 
deviations from the expected quadratic behavior even at rather small Re. 


11 


does not show 
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VIII. CONCLUSIONS 


This paper develops an immersed boundary method that enforces strict rigidity of immersed 
bodies at both zero and finite Reynolds numbers. Unlike existing approaches, we do not rely on 
penalty or splitting approaches, and we instead directly solve a saddle-point system that couples 
the fluid velocity and pressure to the unknown rigidity forces. We developed a physics-inspired 
approximation to the Schur complement (mobility matrix) A4 of the constrained system, based 
on analytical considerations for a continuum fluid model, and demonstrated that this leads to a 
robust preconditioner so long as the immersed boundary markers are kept sufficiently far to ensure 
a well-conditioned mobility matrix. Contrary to common practice, we found that the markers 
should be kept approximately two fluid grid cells apart in rigid-body models in order to obtain 
accurate and stable pointwise estimates for the traction. We tested our method on a number 
of standard test problems in both two and three spatial dimensions, and at both zero and finite 
Reynolds number, and we observed good agreement with theory and literature values. Although 
in this work we focused on rigid bodies, our method can directly be applied to study fluid flow 
around bodies with specified kinematics. For example, it can be used to model the flow around a 
swimming body deforming with a specified gait. We have implemented the method described here 
in the open-source IBAMR software infrastructure [l2] in the hope it will be useful to other users 
of the IB method. 

Another challenge that we did not explore here is the efficient computation of the action of A4 
when there are many markers present; there are many approximate solvers and emerging fast solvers 


we plan to explore in the future. Of course, using dense linear algebra to solve (21) is likely to 


be suboptimal, as these solves have 0{N‘^) memory complexity and 0{N^) time complexity. The 


problem of solving a linear systems similar in structure to (21) appears in many other methods 


for hydrodynamics of suspensions, including Brownian [HIES] and Stokesian m dynamics, the 
method of regularized Stokeslets muHi, computations based on bead models of rigid bodies HZl- 
20], and first-kind boundary integral formulations of Stokes flow m Similar matrices appear 
in static Poisson problems such as electrostatics or reaction-diffusion models [50], and there is a 
substantial ongoing work that can be applied to our problem. Notably, the approximate mobility 
matrix M is dense but has a well-understood low-rank structure that can be exploited. Specifically, 
matrix-vector products A4A can be performed in almost linear time using the Fast Multipole 
Method |58j . If the condition number of A4 is not too large, one can solve linear systems involving 
A4 efficiently using an unpreconditioned Krylov solver. For poorly conditioned cases, however. 
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a good preconditioner based on an approximate factorization of M. is required. In recent years, 
several approximate low-rank factorizations of matrices of this type have been developed [77HZ9], 
and can be used as preconditioners in Krylov methods. We have had reasonable success using a fast 
hierarchically off-diagonal low-rank (HODLR) factorization code developed by Ambikasaran and 
Darve with significant improvement offered by a recently-developed boundary distance low- 
rank approximation m- Preliminary results indicate great promise for the inverse fast multipole 


(iFMM) method [80]; we have been able to use iFMM to solve the system (21) for as many as 
5 • 10^ markers to a relative tolerance of 10“®. These methods are, however, still under active 
development, and a signihcant amount of investigation is necessary to integrate them into the 


method described here. Notably, we only require an approximate solver for (21) and the impact 


of the innacuracy in solving (21) on the overall convergence of the outer Krylov solver needs to be 
assessed. 

The type of linear system we solve here is closely connected to those appearing in implicit 
immersed boundary methods [8TH83] . It is in fact possible to recast the saddle-point problem 
we consider here into a form closely-related to that appearing in implicit IB methods; the Schur 
complement for this system is in Eulerian rather than Lagrangian variables as it was for this work, 
and involves the matrix 


Cy + KS{JS)-^J, 


(28) 


for some constant k that does not need to go to infinity. It may be that geometric multigrid 
methods [S5| developed for implicit IB methods can be applied to the Eulerian Schur complement 


(28). At the same time, techniques developed herein may be useful in the development of more 


efficient implicit IB methods for nearly-rigid bodies. 

Our work is only the first step toward the ultimate goal of developing methods able to handle 
large numbers of rigid bodies in flow. Several computational challenges need to be tackled to 
realize this goal. Firstly, and most importantly, it is crucial to develop a preconditioner for the 
enlarged linear system (??) that appears in the context of freely-moving rigid bodies. An additional 
Schur complement appears when solving this saddle-point problem, and the challenge for future 
work is approximating the body mobility matrix A/" = (/C*A4^^/C) ^. Initial investigations have 
shown great promise in block-diagonal preconditioners with one block per body. In this approach, 
we neglect the hydrodynamic interactions between bodies, but use the mobility approximation 
developed in this work together with dense linear algebra for each body. 


In the marker-based method described in this work, one must adjust the marker spacing to 
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be “neither too small nor too large”. The sensitivity of the solver performance and the numerical 
results to the exact spacing of the markers, which comes from the ill-conditioning of the mobility 
matrix, is one of the key deficiencies of the marker-based representation inherent to the traditional 
IB method. Recently, Griffith and Luo have proposed an alternative IB approach that models 
the deformations and stresses of immersed elastic body using a finite element (FE) representation 
|46| . In their IB/FE approach, the degrees of freedom associated with A are represented on an EE 
mesh that may be coarser than the fluid grid, and the interaction between the fluid grid and body 
mesh is handled by placing IB markers at the numerical quadrature points of the EE mesh. When 
such an approach is generalized to rigid bodies, the conditioning of the mobility becomes much less 
sensitive to the marker spacing. Using a finite-element basis to represent the unknown fluid-body 
interaction force amounts to applying a filter to the marker-based mobility matrix, which is a 
well-known and robust technique to regularize ill-conditioned systems. Specifically, in the context 
of the IB/EE approach, the mobility operator becomes 

Mfe = {JC-^S) 

where is a matrix that contains quadrature weights as well as geometric information about the 
relation between the nodes and quadrature points of the EE mesh. The FE mobility matix A^fe 
is still symmetric, but now can be much smaller because the number of unknowns is equal to the 
number of FE degrees of freedom rather than the number of markers. Even if markers are closely 
spaced, the filtering of the high-frequency modes performed by representing forces in a smooth FE 
basis makes the mobility much better conditioned than for marker-based schemes. Furthermore, 
the mobility matrix, or approximations of it used for preconditioning, will be smaller and thus 
easier to fit in memory. We also expect the resulting method to be more accurate because the 
tractions are represented in a smoother basis. We will explore this promising extension of our 
rigid-body IB methods in future work. 
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Appendix 

Appendix A: Approximating the mobility in three dimensions 


In this appendix we give the details of our empirical hts for the approximations to the fnnctions 


fjsir) and g^(r) in (12) in three spatial dimensions, following the physics-based constraints discussed 
in Section [Tv] To maximize the quality of the fit, we perform separate fits for /? —)> oo (steady Stokes 
flow) and finite /?. We also make an effort to make the fits change smoothly as j3 grows towards 
infinity. 


a. Steady Stokes flow 

Because our numerical computations are done in a periodic domain of length I rather than an 
nnbounded domain, we need to apply a well-known correction to the Oseen tensor IMIETI, 

/oo(h <C r <C 0 ~ (Svrr/r)”^ - 2.84 {67rr]l)~^ . 

From the numerical data, we calculated the normalized functions 

f{x) = {S-n-rjr) (^foo{r) + 2.84/ (bTr??^"^) , (Al) 

g{x) = (87rryr)fi(oo(r), 

where x = r/h is the normalized distance between the markers. As explained previously, we know 
that / ~ (Svrr^r) / (bvrr/a) = 4r/(3a) for a: <C 1 (in practice, markers are never too close to each 
other so we only need the self-mobility, i.e., a; = 0), and that g grows at least quadratically for 
small X (since ^(0) = 0). We also know that / ~ 1 and 5 ~ 1 for large r ^ h. The numerical data 
for the normalized functions f{x) and g{x) are shown in Fig. along with fits to the following 
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semi-empirical rational functions, 


fix) = < 

gix) = 


{3a) / {4:h)+box'^ 


bixe 


—b 2 X 


+ 


bzx^+x'^ 

l+b^x'^-l-x'^ 


if X < 0.8, 
if X > 0.8, 


(A 2 ) 


65 -I- 60x2 -I- x^ ■ 

As the figure shows, the numerical data are well described by these formulas, and there is only 
small scatter of the numerical data around the fit, indicating approximate discrete translational 


and rotational invariance We also obtain a reasonable agreement with the RPY tensor (16) 
approximation; however, as expected, the empirical fits yield a better match to the data. 


b. Nonzero Reynolds numbers 


For hnite /3, we consider separately the case r = 0 (giving the diagonal elements M-a) and 


T > O.l/i (giving the off-diagonal elements). For r 


to (18), 


0 we use an empirical fit designed to conform 


^ yh/isjo) ^ 1 -F Zi^/P -I- Z2I3 

P zo + zsP -I- 6Tr{a/h)z2P‘^ ’ 

5 / 3 ( 0 ) = 0 , 


(A3) 


where zi — Z 3 are coefficients obtained by fitting the numerical data for the self mobility for different 
p. Note that zq = 2hP j {SVm) is fixed by the inviscid condition Also note that as /? —>■ 00 , 

our fit obeys the correct Stokes limit, 


We show the empirical fit for (pQ (P) in Fig. 16 in Appendix [B| 

For nonzero r, we introduce normalized functions fp and via 


P 




rjh Attx^ 


9pi'^) = in ■ 


gh dvrx^ 


fyix), 

gyix), 


(A4) 


11 


Most of the scatter comes from the finite size of the periodic box and can be explained using a known periodic 
correction to the RPY tensor | 84| . 
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where x = r/h is the normalized distance, and fi/rj = /S.t/ph?. For finite /?, we know that fpix ^ 


y/P) ~ 5 , 3 ( 2 ; 3> \/~P) « 1 according to (20). As /3 —>■ 00 , we want to reach the Stokes limit 

foo (X > 1) ^ , 

x2 

Qoo (x > 1 ) , 


(A5) 


and for finite /3, we want the viscous contribution to decay as exp (—x/ for some constant 

C that should be close to unity. Furthermore, we would like to ensure continuity near the origin 
with the fit for r = 0 , 

fpix 0 ) -A-KX^ifo (/?). 

A fitting formula that obeys these conditions that we find to work well for r > O.lh is 


—dvrx^ + 04 x^ — x'^e^ 

fp{ ) ^oiP) 1 + a4X®(/Jo(/3) 


+ 


+ 


05x^6 + ajx^ 

1 + osx^ + agx^ ’ 


(A6) 


X® + x'^e^ boxly/Jj) 


gp{x) ^ ^^^2 _j_ 5^2.3 _|_ 5^3,4 _|_ ’ 


where ag-ag and 60-^5 are empirical coefficients. It is important to emphasize that (A 6 ) was chosen 


in large part based on empirical trial and error. Many other alternatives exist. For example. 


one could use the analytical Brinkmanlet (15) for sufficiently large distances and then add short- 


ranged corrections for nearby markers. Alternatively, one could first subtract the inviscid part 
/o(r) and go{r) and then fit the viscous contribution only. As discussed above, ideally the fits 
would be constrained to guarantee an SPD approximate mobility matrix, but this seems difficult 
to accomplish in practice. 


We computed the fitting coefficients in (A 6 ) for /3 G {0, 0.1, 0.25, 0.5, 1, 10, 100, 1000}; the 
coefficients for other values in the range 0 < /3 < 1000 are interpolated using linear interpolation, 
and f3 > 1000 is treated using the steady Stokes fitting. We see a good match between the 
numerical data and our empirical fits in Fig. pAj with good translational and rotational invariance 
(i.e., relatively small scatter of the numerical points around the fits). 


Appendix B: Approximating the mobility in two dimensions 


To construct empirical approximations to the functions fp^r) and 5 ^ 3 (r) in (12) in two spatial 
dimensions, we follow the same approach as we did for three dimensions in Appendixj^ Specifically, 
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Figure 14: Normalized fitting functions fp{x) (left) and gp{x) (right) at finite /3 in three dimensions for the 
6-point kernel, for different values of the viscous CFL number (see legend). Symbols are numerical data 


obtained by using a 256^ periodic fluid grid, and dashed lines show the best fit of the form (A6). 


we first discuss the known asymptotic behavior of these functions at short and large distances, and 
use this to guide the construction of empirical fitting formulas. 


1. Physical constraints 


In two dimensions, we need to modify (18) to agree with (17) for small /3. For d = 2, Vm = Cyh'^ 
and /o(0) ^ 13/rj so that we use the fit 


U (0) = and 5/3(0) = 0, 


7] 


(Bl) 


where C'(/?) has the same asymptotic scaling as in three dimensions and is obtained from empirical 


fits (see Fig. 16). A key difference exists between two and three spatial dimensions in the limit 
Re —)■ 0. For steady Stokes flow in a square two dimensional periodic domain, the Green’s function 
diverges logarithmically with the system size 1. Therefore, it is not possible to write a formula for 
the asymptotic behavior at large distances for an infinite system. Instead, we must subtract the 
divergent piece to get a well-defined answer. The standard Green’s function for Stokes flow in two 


dimensions has logarithmic growth at infinity, which suggests that (|19|) should be replaced by 

1 


fooir > /i) - foo{r = 0) Ri and goo{r > /i) ~ -—. 

dvrr/ 4:7rr] 


For inviscid flow we should replace (20) by the field of a dipole in two dimensions. 

At 


fo{r > /i) Ri - 


At 

27rpr"^ 


and go{r ^ h) 


7rpr 


2 * 


(B2) 


(B3) 


In two dimensions the solutions of the Brinkmann equation (14) are analytically complicated and 


involve special functions. Even without solving these equations, however, physical scaling suggests 
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that the same physical length scale h^/]3 should enter, in particular, the viscous corrections should 
decay to zero exponentially fast with /ly^. 


2. Empirical fits 

We have used the analytical results above to construct empirical fitting formulas that have the 
correct asymptotic behavior, as we now explain in more detail. 


a. Steady Stokes flow 


In two dimensions, steady Stokes flow (/3 —oo) is not well behaved because the Green’s function 
does not decay sufficiently rapidly (Stokes paradox). This makes the mobility an essentially dense 
matrix that is sensitive to boundary conditions and difficnlt to approximate. Nevertheless, we have 


used a periodic system to fit empirical data based on the theory (B2). The diagonal value /oo(0) 
diverges logarithmically with the system size L for periodic boundaries. Specifically, for a square 
unit cell of length I ^ h, it is known that m 


and this relation defines the effective hydrodynamic radius of a marker a (note that a/h is a 
universal value for a given spatial discretization, as it is in three dimensions). Since the precise 
form depends on boundary conditions and is not known in general, we treat /oo(0) as an input 
parameter. 

From the numerical data, we calculated the normalized functions 


/(x) = - (47rr/) (/oo(x) -/oo(0)) (B4) 

g{x) = (47rr/)5oo(x), 


where x = r/h is the normalized distance between the markers. Observe that from (B2) we know 
that /(x S> 1) ~ Inx and g{x ;:?> 1) ~ 1. For the normalized functions, we use the fits 


fix) 

six) 


oox^ + aix^ + a2X^ In x 
1 + oax + 04x2 + a2X^ ’ 
6ox^ + bix^ 

1 + 62X + 63x2 + bix ^ 


(B5) 


Numerical results and empirical fits for /(x) and g{x) are shown in Fig. 


15 


While the numerical 


data do conform to the theoretical asymptotic behavior, there is snbstantial scatter for larger 


distances because of the strong sensitivity to the boundaries. 
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• L=256 
L=512 


Figure 15: Empirical fits (lines) to numerical data (symbols) for /(x) (left) and g{x) (right), for the 6-point 
kernel in two dimensions, obtained using a periodic system of either 256^ or 512^ grid cells. Observe that 
both follow the correct asymptotic behavior at large distances, with scatter dominated by boundary effects. 


b. Nonzero Reynolds numbers 


For r = 0, we use a fitting formula in agreement with (Bl), 


‘fo (/ 3 ) = 


_ zo + zi/3HogiP) 
f3 1 -b Z2(3 -b Z3/32 -b 2:4/34 ’ 
9/3(0) = 0 , 


(B6) 


where 2:1 — 2:4 are coefficients (obtained by fitting for each kernel data over a range of /3’s) and zq 
is fixed from the inviscid condition The empirical fit for tpo (/3) is shown in Fig. Note 

that for finite /3, one must ensure that the system size used to tabulate the values of fjs and gp is 
sufficiently large, I hy/]3. 



Figure 16: Empirical fit for ipo (/3) as a function of /3 for different vaues of f3. (Left) Three dimensions, 256^ 
grid. (Right) Two dimensions, 512^ grid. 
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For r > O.lh we introduce normalized functions and gjs via 


fl3{r) = - 
9^ir) = 




27vr]x^ 




(B7) 


13 


TTrjx^ 


■9l3{3 


where in the inviscid case we take /3/?? = /S.t/ph?, and x = r/h is the normalized distance. For 


finite /3, we know that fp{x >> ^/]3) ~ S> ^/]3) ^ 1 according to (B3). The numerical data is 
fitted with the empirical fitting functions 


~ , x^ ln(x) -iiL 

fp{x) = — - ^ + 


9fiix) = 


/3(ao + 2x) 

^3 


_ P2" 

e + 


+ a2X^ + asx^ 

1 + + b2X^ + asx'^ ’ 

^3 


(B8) 


f5{cQ + Ax) ' e P3®(ci + C 2 X + caa:^) + ’ 

as shown in Fig. dzi Here 09 — 03, pi — P 3 , 61 — 63 , cq — C 3 are empirical coefficients, computed 
by fitting numerical dats for /3 in {0, 0.1, 0.25, 0.5, 1.0, 5.0, 10.0}. Intermediate values in the 
range 0 < (3 < 10.0 are interpolated using linear interpolation, and larger /3’s are handled using 


the steady Stokes fit (B5). 



-0.5 


-p=io 


-p=io 


Figure 17: Empirical fitting of fp{r) and gp^r) in two dimensions for different values of j5 for the 6-point 
kernel, 512^ grid. 


Appendix C: Stokes flow between two concentric spheres 

Consider steady Stokes flow around a rigid shell or sphere of radius a, placed in a centered 
position inside another spherical shell or cavity of radius 6 = a/A. We consider the case when 
the outer shell is moving with velocity V and the inner shell is at rest, and for simplicity set the 
viscosity to unity, p = 1. Brenner |85j gives the drag force on the inner sphere for no slip boundary 
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conditions as 


F = —GnaVK, 


(Cl) 


where 


1 --^^ 1 . 9 , 5,3 9,5 ,6 

K —- cnid. ex — 1 — —A-|-—A — —A “hA. 

a 4 2 4 


Let us denote the constants 


15C A^-A^ 


B = 


D = 


4a2 

a ’ 

2)V a 1 

-A5 

2 

a 

V 1 + 

5 x3 9\5 

2 

a 

Va^ 1 

-A3 

4 

a 


The velocity in the region between the two spherical shells can be obtained from the expressions 
given by Brenner as 

Vr = — COS 6 * (—r^ — B- + 2C + 211^ ) ; 

V 5 r r'^ / 


, /yl n Bl _ 1 , 

U0 = sin0( -r --- + 2 C'-Z 1 ^ ) , 


and the pressure is 


= 0 , 


^cosH ^ . 

TT = TToo + f^B -2- Z^Ar cost), 

where ttoo = 0 since we impose that the pressure have mean zero to remove the null mode for 
pressure. In spherical coordinates, with the symmetry axes aligned with the direction of the flow, 
the traction on the surface of the inner sphere, which is the jump in the stress across the inner 
shell, is 


X = a ■ n 


H cos 9 




f + /i sin 0 



e, 


where r = (sin 9 cos (f>, sin 9 sin (f), cos 9) and 6 = (cos 9 cos (f), cos 9 sin 0, — sin 9). 
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Appendix D: Imposing Physical Boundary Conditions 

The local averaging and spreading operators have to be modified near physical boundaries, 
specifically, when the support of the kernel Sa overlaps with a boundary. A proposal for how to do 
that has been developed by Yeo and Maxey [86] , and an alternative proposal has been developed in 
the context of the immersed boundary method by Griffith et al. 1HZ|. Here we have chosen to use 
the former approach because of its simplicity and the fact that it is independent of the kernel, as 
well as the fact that it ensures that the interpolated velocity strictly vanishes at a no-slip boundary; 
this ensures that the mobility of a marker is a monotonically decreasing function as it approaches a 
no-slip boundary. Since the description in |86| is limited to steady Stokes flow and a single no-slip 
boundary, we give here an algebraic formulation that extends to a variety of boundary conditions; 
this formulation is implemented in the IBAMR library and used in the examples in this paper in 
non-periodic domains. 

The basic idea in the handling spreading and interpolation near boundaries is to use the standard 
IB kernel functions in a domain extended with sufficiently many ghost cells so that the support of 
all kernels is strictly within the extended domain. For interpolation, we first fill ghost cells and 
then interpolate as usual using the ghost cell values. For spreading, we take the adjoint operator, 
which basically means that we first spread to the extended domain including ghost cells in the 
usual manner, and then we accumulate the value spread to the ghost cell in the corresponding 
interior grid point, using the same weight (coefficient) that was used when filling ghost cells for 
the purposes of interpolation. 

This process requires a consistent method for filling ghost cells, that is, for extending a (cell- 
centered or staggered) field u from the interior to the extended domain. In general, this will be an 
affine linear mapping of the form 

"^ext — C, 

where E is an extension matrix and c encodes inhomogeneous boundary conditions. Let us denote 
with jT’o the standard IB interpolation operator that interpolates an extended field at a position 
inside the interior of the domain. The interpolated value in the presence of physical boundary 
conditions is then given by the affine linear mapping 

BC (’^int) — O^ext — 0-^’^int T 0^- 

The corresponding spreading operation is defined to be the adjoint of JF bc for homogeneous bound¬ 
ary conditions, as this ensures energy conservation in the absence of boundary forcing. Specifically, 
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we use 


Sbc = = E^So. 

The specific form of the extension operator E used in our implementation is based on linear 
extrapolation to a given ghost point based on the corresponding value in the interior and the values 
at the boundary as specified in the boundary conditions. Specifically, for homogeneous Neumann 
conditions we do a mirror image Ughost = —Uint, while for homogeneous Dirichlet boundary con¬ 
ditions, such as no slip boundaries, we simply do a mirror inversion Ughost = —Uint- This makes 
our implementation exactly identical to that proposed by Yeo and Maxey [86] in the context of 
the FCM method. One can think of this approach to no-slip boundaries as taking an inverted 
mirror image of the portion of the kernel outside of the domain |86| . Note that the same E is used 
to implement boundary conditions both in the fluid solver and when interpolating/spreading near 
boundaries; this greatly simplifies the implementation without lowering the second-order accuracy 
of the fluid solver |32| . In our implementation, we use simple transpose for spreading which is the 
adjoint E* with respect to the standard inner product. 
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